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Abstract 

A description of the IntraNuclear Cascade (INC), preequilibrium, evaporation, fission, 
coalescence, and Fermi breakup models used by the latest versions of our CEM03.03 and 
LAQGSM03.03 event generators is presented, with a focus on our most recent developments of 
these models. The recently developed "S" and "G" versions of our codes, that consider multi- 
fragmentation of nuclei formed after the preequilibrium stage of reactions when their excitation 
energy is above 2A MeV using the Statistical Multifragmentation Model (SMM) code by Botv- 
ina et al. ("S" stands for SMM) and the fission-like binary-decay model GEMINI by Charity 
("G" stands for GEMINI), respectively, are briefiy described as well. Examples of benchmark- 
ing our models against a large variety of experimental data on particle-particle, particle-nucleus, 
and nucleus-nucleus reactions are presented. Open questions on reaction mechanisms and future 
necessary work are outlined. 
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1. Introduction 



Following an increased interest in intermediate- and high-energy nuclear data in relation 
to such projects as the Accelerator Transmutation of nuclear Wastes (ATW), Accelerator Pro- 
duction of Tritium (APT), Spallation Neutron Source (SNS), Rare Isotope Accelerator (RIA), 
Proton Radiography (PRAD) as a radiographic probe for the Advanced Hydro-test Facility, 
NASA needs, and others, the US Department of Energy has supported during the last decade 
our work on the development of improved versions of the Cascade-Exciton Model (CEM) and 
of the Los Alamos version of the Quark Gluon String Model (LAQGSM) which has led to our 
intermediate- and high-energy event generators CEM03.03 and LAQGSM03.03, respectively, 
and their modifications described below. 

The main focus of our workshop is nucleon-induced reactions for Accelerator Driven Systems 
(ADS) and spallation neutron sources (up to 2-3 GeV), usually described well enough by our 
intermediate-energy code CEM03.03 without a need to use our high-energy code LAQGSM03.03. 
This is why we focus below mostly on CEM03.03 and its modifications. However, as discussed 
below in Section 3.1, CEM does not consider the so-called "trawling" effect (depletion of target 
nucleons during a cascade), therefore does not describe well reactions on very light nuclei like 
C at incident energies above about 1 GeV. Therefore, in transport codes that use both CEM 
and our high-energy code LAQGSM as event generators, we recommend simulating nuclear 
reactions with CEM at incident energies up to about 1 GeV for light nuclei like C and up to 
about 5 GeV for actinide nuclei, and to switch to simulations by LAQGSM, which does consider 
the "trawling" effect, at higher energies of transported particles. This is the reason we have 
included in the present lectures a brief description of LAQGSM as well. Even at energies of 
ADS applications below about 3 GeV, we recommend using LAQGSM instead of CEM in the 
case of light target nuclei. 

The Cascade-Exciton Model (CEM) of nuclear reactions was proposed almost 30 years ago 
at the Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 
Dubna, USSR by Gudima, Mashnik, and Toneev [Il|2]. It is based on the standard (non time- 
dependent) Dubna IntraNuclear Cascade (INC) P g] and the Modified Exciton Model (MEM) 
[3, E] • CEM was extended later to consider photonuclear reactions [7J and to describe fission 
cross sections using different options for nuclear masses, fission barriers, and level densities [8] 
and its 1995 version, CEM95, was released to the public via NEA/OECD, Paris as the code 
IAEA1247, and via the Radiation Safety Information Computational Center (RSICC) at Oak 
Ridge, USA, as the RSICC code package PSR-357 [S]. 

The International Code Comparison for Intermediate Energy Nuclear Data pUl [TT] orga- 
nized during 1993-1994 at NEA/OECD in Paris to address the subject of codes and models 
used to calculate nuclear reactions from 20 to 1600 MeV showed that CEM95 (more exactly, 
CEM92m, which is almost identical in its physics content to the publicly released CEM95 ver- 
sion) had one of the best predictive powers to describe nucleon-induced reactions at energies 
above about 150 MeV when compared to other models and codes available at that time. 

The code LAQGSM03.03 described briefly below is the latest modification ^ of LAQGSM [13], 
which in its turn is an improvement of the Quark-Gluon String Model (QGSM) [T3]. It describes 
reactions induced by both particles and nuclei at incident energies up to about 1 TeV/nucleon, 
generally, as a three-stage process: IntraNuclear Cascade (INC), followed by preequilibrium 
emission of particles during the equilibration of the excited residual nuclei formed after the 
INC, followed by evaporation of particles from and/or fission of the compound nuclei. 
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CEM95 and/or its predecessors and successors like CEM97 [I5l[l6], CEM2k [17], 
CEM2k+GEM2 [l8]-|io], CEM03 |2ll[22], including the latest version available to the pubhc 
from RSICC and NEA OECD, CEM03.01 |23|, as well as different versions of LAQGSM and 
its predecessor QGSM are used as stand-alone codes to study different nuclear reactions for 
applications and fundamental nuclear physics (see, e.g., and references therein). Parts 

of different versions of the GEM, LAQGSM, or QGSM codes are used in many other stand- 
alone codes, like PICA95 |13], PICAS CASCADO [35], CAMO g6], MCFX |17j, 
ECM t48j, and NUCLEUS |19]. Different versions of GEM and LAQGSM, or of the older 
QGSM, are incorporated wholly, or in part in many transport codes used in a variety of appli- 
cations, like CASCADE [SU], MARS [51], MCNPX [S2] GEANT4 [SSIIM], SHIELD [52], 
RTS&T [56], SONET [57], CALOR [58], HETC-3STEP p], CASCADE/INPE [60], 
HADRON [61], CASCADE-2004 [62], and others. 

In these lectures, we present a brief description of all models, approximations, and system- 
atics used in the latest versions of our GEM and LAQGSM event generators. 

2. A Brief Survey of CEM and LAQGSM Physics 

The basic version of both our GEM and LAQGSM event generators is the so-called "03.01" 
version, namely GEM03.01 [23] and LAQGSM03.01 [21]. The GEM03.01 code calculates nuclear 
reactions induced by nucleons, pions, and photons. It assumes that the reactions occur generally 
in three stages. The first stage is the IntraNuclear Gascade (ING), in which primary particles 
can be re-scattered and produce secondary particles several times prior to absorption by, or 
escape from the nucleus. When the cascade stage of a reaction is completed, GEM03.01 uses 
the coalescence model to "create" high-energy d, t, ^He, and ^He by final-state interactions 
among emitted cascade nucleons, already outside of the target. The emission of the cascade 
particles determines the particle-hole configuration, Z, A, and the excitation energy that is the 
starting point for the second, preequilibrium stage of the reaction. The subsequent relaxation 
of the nuclear excitation is treated in terms of an improved version of the modified exciton 
model of preequilibrium decay followed by the equilibrium evaporation/fission stage of the 
reaction. Generally, all four components may contribute to experimentally measured particle 
spectra and other distributions. But if the residual nuclei after the ING have atomic numbers 
with A < 12, GEM03.01 uses the Fermi breakup model to calculate their further disintegration 
instead of using the preequilibrium and evaporation models. Fermi breakup is much faster 
to calculate and gives results very similar to the continuation of the more detailed models to 
much lighter nuclei. As already mentioned in the Introduction, LAQGSM03.01 also describes 
nuclear reactions, generally, as a three-stage process: IntraNuclear Gascade (ING), followed by 
preequilibrium emission of particles during the equilibration of the excited residual nuclei formed 
after the ING, followed by evaporation of particles from or fission of the compound nuclei. 
LAQGSM was developed with a primary focus on describing reactions induced by nuclei, as 
well as induced by most elementary particles, at high energies, up to about 1 TeV/nucleon. The 
ING of LAQGSM is completely different from the one in GEM. LAQGSM03.01 also considers 
Fermi breakup of nuclei with A < 12 produced after the cascade, and the coalescence model to 
"produce" high-energy d, t, ^He, and ^He from nucleons emitted during the ING. 

The main difference of the following, so-called "03.02" versions of GEM and LAQGSM from 
the basic "03.01" versions is that the latter use the Fermi breakup model to calculate the 
disintegration of light nuclei instead of using the preequilibrium and evaporation models only 
after the ING, when the excited nuclei after the ING have a mass number A < 12, but do not 
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use the Fermi breakup model at the preequihbrium, evaporation, and fission stages, when, due 
to emission of preequihbrium particles or due to evaporation or to a very asymmetric fission, 
we get an excited nucleus or a fission fragment with A < 12. This problem was solved in the 
03.02 versions of CEM and LAQGSM, where the Fermi breakup model is used at any stage of 
a reaction, when we get an excited nucleus with A < 12. 

In addition, the routines that describe the Fermi breakup model in the basic 03.01 version 
of our codes were written some twenty years ago in the group of Prof. Barashenkov at JINR, 
Dubna, Russia, and are far from being perfect, though they are quite reliable and are still 
used currently without any changes in some transport codes. First, these routines allow in 
rare cases production of some light unstable fragments like ^He, ^Li, ^Be, ^B, etc., as a result 
of a breakup of some light excited nuclei. Second, these routines allowed in some very rare 
cases even production of "neutron stars" (or "proton stars"), i.e., residual "nuclei" produced 
via Fermi breakup that consist of only neutrons (or only protons). Lastly, in some very rare 
cases, these routines could even crash the code, due to cases of 0/0. All these problems of the 
Fermi breakup model routines are addressed and solved in the 03.02 version of our codes |63j . 
Several bugs are also fixed in 03.02 in comparision with its predecessor. On the whole, the 

03.02 versions describe nuclear reactions on intermediate and light nuclei, and production of 
fragments heavier than ^He from heavy targets much better than their predecessors, almost do 
not produce any unstable unphysical final products, and are free of the fixed bugs. 

However, even after solving these problems and after implementing the improved Fermi 
breakup model into CEM03.02 and LAQGSM03.02 153], in some very rare cases, our event 
generators still could produce some unstable products via very asymmetric fission, when the 
excitation energy of such fragments was below 3 MeV and they were not checked and not 
disintegrated with the Fermi breakup model (see details in [12]). This problem was addressed 
in the 03.03 versions of our codes, where we force such unstable products to disintegrate via 
Fermi breakup independently of their excitation energy. Several more bugs were fixed on the 

03.03 version as well. A schematic outline of a nuclear reaction calculation by CEM03.03 or 
LAQGSM03.03 is shown in Fig. 1. We emphasize that the occurrence of these problems even 
in the 03.01 versions is quite rare, allowing stand-alone calculations of many nuclear reactions 
to proceed without problems, but are unacceptable when the event generators are used inside 
transport codes doing large-scale simulations. 

In the following Sections, we highlight the main assumptions of the models contained in 
CEM and LAQGSM. 

3. The Intranuclear Cascade Mechanism 

The inelastic interaction of a high-energy particle with a nucleus, and even more the colli- 
sions of two nuclei, is a very complex and multi-faceted phenomenon whose analytical descrip- 
tion encounters considerable difficulties [3l H] . In recent years calculations of such interactions 
have been carried out by statistical simulations using the Monte-Carlo method. 

The INC approach was apparently first developed by Goldberger [Blj, who in turn based 
his work on the ideas of Heisenberg and Serber [55], who regarded intranuclear cascades as a 
series of successive quasi-free collisions of the fast primary particle with the individual nucleons 
of the nucleus. 

Let us recall here the main basic assumptions of the INC, following [31 H]. The main condition 
for the applicability of the intranuclear-cascade model is that the DeBroglie wavelength A of 
the particles participating in the interaction be sufficiently small: It is necessary that for 
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Figure 1: Flow chart of nuclear-reaction calculations by CEM03.03 and LAQGSM03.03. 



most of these particles A be less than the average distance between the intranuclear nucleons 
A ~ 10~^^ cm. Only in this case does the particle acquire quasi-classical features and can we 
speak approximately of particle trajectory and two-particle collisions inside the nucleus. It is 
clear that for this to be the case the primary particle kinetic energy T must be greater than 
several tens of MeV. 

Another important condition for applicability of the INC is the requirement that the time 
in which an individual two-particle intranuclear collision occurs on the average, r ~ 10"^'^ sec, 
be less than the time interval between two such consecutive interactions 

At = l/c > AnRySAac > 3 ■ 10'^^ a (mb) sec, 

where / is the mean range of the cascade particle before the interaction, c is the velocity of 
light, R = tqA^^^ is the mean radius of the nucleus, and a is the cross section for interaction 
with an intranuclear nucleon. This permits the interaction of the incident particle with the 
nucleus to be reduced to a set of individual statistically independent intranuclear collisions. 

The requirement r < At is equivalent to the requirement that the intranuclear interaction 
cross section be sufficient small: a < 100.^ mb, where the coefficient ~ 1. 



5 



Since the energy of the particles participating in the cascade is rather large — as a rule signif- 
icantly greater than the binding energy of the intranuclear nucleons — the same characteristics 
can be used for interaction of cascade particles inside the nucleus as for the interaction of free 
particles. The effect of other intranuclear nucleons is taken into account by introduction of 
some average potential V, and also by the action of the Pauli principle.^ 

We can say that a high-energy particle which has entered the nucleus passes through a gas 
of free nucleons, producing a cascade (avalanche) of secondary particles. A fraction of these 
secondary particles leaves the nucleus, and the remaining fraction is absorbed, exciting the 
nucleus to some energy E*. 

Following [21 S], after the choice of a nuclear model and an algorithm for determination of 
the elementary particles involved in the INC with the intranuclear nucleons (for this purpose 
it is necessary to store in the computer memory also the values of the integrated cross sections 
for elastic and inelastic interactions o"ez(T) and crj„(T)), calculation of the intranuclear cascade 
is carried out according to the scheme shown in Fig. 2. The turquoise boxes 1, 2, 4, 5, 8-10, 
12, and 14 in the diagram denote operations which are definite logically closed parts of the 
INC program. The yellow boxes 3, 6, 7, 11, and 13 denote logical operations which control the 
various branchings of the program (transfer conditions). 

Box 1 takes into account the change in primary-particle momentum due to the effect of the 
intranuclear potential and to refraction and reflection of the DeBroglie wave of the particle at 
the nuclear boundary. 

In box 2 is chosen the momentum and isospin (proton or neutron) of the intranuclear nucleon 
with which the interaction occurs (for brevity we will call this nucleon the partner), and from 
the given elementary cross section atot{t) = Cez(i) + o"m(t) (where t is the relative energy of 
the primary particle and the partner taking part in the intranuclear motion) the mean free 
path of the particle in nuclear matter L = L{atot) is calculated and the point of interaction is 
determined. 

Box 3 tests whether this point of interaction is inside the nucleus. If it is not, then the 
particle is assumed to have passed through the nucleus without interaction. The ratio of the 
number of such particles to the total number of interactions considered with the nucleus Ntot 
obviously characterizes the reaction cross section ain{t). 

If the point of interaction is inside the nucleus, then the type of interaction: elastic or 
inelastic, is determined from the known cross sections aei{t) and crj„(t) in box 4. 

In box 5 the secondary-particle characteristics are determined in accordance with the type 
of interaction selected (the nature, number, energy, and the emission angle). 

Box 6 is a test of whether the Pauli principle is satisfied. Interactions which do not satisfy 
this principle are considered forbidden and the particle trajectory is followed beyond the point 
of the forbidden interaction. 

In box 7 the particle energy T is compared with some initially specified cutoff energy Tcut 
which determinates whether this particle is sufficiently energetic (T > Tcut) to take further part 
in development of the intranuclear cascade or whether its energy is so small (T < T^ut that the 
particle is simply absorbed by the nucleus. In the first case the particle is followed further as 
was described above. (For this the parameters of all cascade particles with energy T > T^-ut cire 
stored in the memory in box 8 and later the cascade calculation is repeated for each of them 

"'^The nucleus is considered to be a degenerate Fermi gas of nucleons enclosed in the nuclear volume. According 
to the Pauli principle the nucleons, after an intranuclear collision, must have energy above the Fermi energy; 
otherwise such an interaction is forbidden. The action of the Pauli principle leads in effect to an increase of the 
mean free path of fast particles inside the nucleus. 
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Figure 2: Flow chart for the intranuclear cascade calculation. 



in turn by going to boxes 9 and 2.) In the second case the INC treatment of this particle is 
terminated; if this particle is a nucleon, in box 10 it contributes to the energy of the residual 
nucleus and become an exciton to be further treated by the Modified-Exciton Model in box 12. 

The calculation is carried out until all particles are absorbed or leave the nucleus. The 
operations in boxes 8, 9, and 11 are responsible for this. If the history of one particle which 
entered the nucleus had been completed {i.e., if the computer memory is empty; see box 11), 
a possible preequilibrium, followed by evaporation/fission, or/and Fermi breakup stage of this 
event is simulated in box 12 until the excitation energy of the residual nucleus is below the 
binding energy of a neutron or other particles that could be emitted from this nucleus, then, 
the history of the next particle {i.e., next "event") is simulated (boxes 13 and 14), and so forth, 
until all events are simulated and we get the needed statistics. 

Any cascade calculation at not very high energies where it is still possible to neglect many- 
particle interactions and the change in density of the intranuclear nucleons can be fitted into the 
general scheme shown in Fig. 2. The specific form of the box operations and their complexity 
are determined by the choice of the nuclear model and by the number and variety of elemen- 
tary processes which it is considered necessary to take into account in a given calculation. The 
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individual boxes can be studied in more detail in Refs. [661 EZ], as well as in the quite old, 
but still one of the best monographs on the INC and other nuclear reaction models we highly 
recommend to readers interested in details of high-energy nuclear reactions [3j . Some specific 
details on the INC of CEM and LAQGSM are provided in the following Sections 3.1 and 3.2. 



3.1. The INC of CEM03.03 



The intranuclear cascade model in CEM03.03 is based on the standard (non-time-dependent) 
version of the Dubna cascade model [21 IH ESI EZj- AH the cascade calculations are carried 
out in a three-dimensional geometry. The nuclear matter density p(r) is described by a Fermi 
distribution with two parameters taken from the analysis of electron-nucleus scattering, namely 

Pi.r) = ppir) + pnir) = po{l + exp[{r - c)/a]} , (1) 

where c = 1.07A^^^ fm, A is the mass number of the target, and a = 0.545 fm. For simplicity, 
the target nucleus is divided by concentric spheres into seven zones in which the nuclear density 
is considered to be constant. The energy spectrum of the target nucleons is estimated in the 
perfect Fermi-gas approximation with the local Fermi energy Tp{r) = /i^[37r^p(r)]^/^/(2mAr), 
where is the nucleon mass. An example of the nucleon density and the Fermi energy used 
by CEM03.01 to calculate nuclear reactions on ^osp^ jg shown in Fig. 3. 

The influence of intranuclear nucleons on the incoming projectile is taken into account by 
adding to its laboratory kinetic energy an effective real potential V, as well as by considering 
the Pauli principle which forbids a number of intranuclear collisions and effectively increases 
the mean free path of cascade particles inside the target. For incident nucleons V = V^lr) = 
Tp{r) + e, where Tp{r) is the corresponding Fermi energy and e is the binding energy of the 
nucleons. For pions, CEM03.01 uses a square-well nuclear potential with the depth — 25 
MeV, independently of the nucleus and pion energy, as was done in the initial Dubna INC OH). 
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Figure 3: Examples of the nucleon density and the Fermi energy used by CEM03.01 to calculate 
nuclear reactions on ^''^Pb. 

The interaction of the incident particle with the nucleus is approximated as a series of 
successive quasi-free collisions of the fast cascade particles (A^, vr, or 7) with intranuclear 
nucleons: 



NN NN, NN ^ TcNN, NN ^ 7^1, 

■nN^nN, tcN tti, ■ ■ ■ ,7riN {i > 2) . 



(2) 
(3) 
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In the case of pions, besides the elementary processes (3), CEM03.01 also takes into account 
pion absorption on nucleon pairs 

ttNN NN. (4) 

The momenta of the two nucleons participating in the absorption are chosen randomly from 
the Fermi distribution, and the pion energy is distributed equally between these nucleons in 
the center-of-mass system of the three particles participating in the absorption. The direction 
of motion of the resultant nucleons in this system is taken as isotropically distributed in space. 
The effective cross section for absorption is related (but not equal) to the experimental cross 
sections for pion absorption by deuterons. 

In the case of photonuclear reactions [22], CEM03.01 follows the ideas of the photonuclear 
version of the Dubna INC proposed initially 35 years ago by one of us (KKG) in collaboration 
with Iljinov and Toneev [68] to describe photonuclear reactions at energies above the Giant 
Dipole Resonance (GDR) region [69j. [At photon energies = 10-40 MeV, the DeBroglie 
wavelength ^ is of the order of 20-5 fm, greater than the average inter-nucleonic distance in 
the nucleus; the photons interact with the nuclear dipole resonance as a whole, thus the INC 
is not applicable.] Below the pion-production threshold, the Dubna INC considers absorption 
of photons on only "quasi-deuteron" pairs according to the Levinger model [70j : 

r ZjA-Z) 

where A and Z are the mass and charge numbers of the nucleus, L ^ 10, and cXyd is the total 
photo-absorption cross section on deuterons as defined from experimental data. 

At photon energies above the pion-production threshold, the Dubna INC considers produc- 
tion of one or two pions; the specific mode of the reaction is chosen by the Monte-Carlo method 
according to the partial cross sections (defined from available experimental data): 

7 + p — > p + 7r° , (6) 
^ n + 7r+ , (7) 

^ J9 + 7r+ + TT" , (8) 

^ J9 + 7r° + 7r° , (9) 
n + 7r+ + 7r° . (10) 

The cross sections of 7 + n interactions are derived from consideration of isotopic invariance, 
i.e. it is assumed that (7(7 + n) = 17(7 + p). The Compton effect on intranuclear nucleons is 
neglected, as its cross section is less than ^ 2% of other reaction modes (see, e.g. Fig. 6.13 in 
Ref. [TT]). The Dubna INC does not consider processes involving production of three and more 
pions; this limits the model's applicability to photon energies < 1.5 GeV [for higher than 
the threshold for three-pion production, the sum of the cross sections (8)-(10) is assumed to be 
equal to the difference between the total inelastic 7 + p cross section and the sum of the cross 
sections of the two-body reactions (6)-(7)]. 

The integral cross sections for the free NN, ttN, and 7iV interactions (2)-(10) are approxi- 
mated in the Dubna INC model [3] used in CEM95 [9| and its predecessors using a special algo- 
rithm of interpolation/extrapolation through a number of picked points, mapping as well as pos- 
sible the experimental data. This was done very accurately by the group of Prof. Barashenkov 
using all experimental data available at that time, about 35 years ago. Currently the exper- 
imental data on cross sections is much more complete than at that time; therefore we have 
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Figure 4: Energy dependence of the np, vr+p, and 7p total cross sections and of the np — ^ ppvr", 
Ti^p — » 7r+p7r°, and 7p mi^TT^ ones. Experimental points are from our compilation |15j . 
Solid hues are results using our new approximations; dashed lines show the standard Dubna 
INC approximations [3] used in CEM95 [9]. 
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Figure 5: Comparison of eight experimental total 7 + p{d) cross sections with the old ap- 
proximations used in the Dubna INC [3] and with our approximations incorporated into the 
CEM03.01 and LAQGSM03.01 codes. The red curve gives the code results using parabolic 
interpolation, while the blue solid curve uses linear interpolation between our tabulated points. 
Where no blue curve is visible, it is coincident with the red curve. References to experimen- 
tal data shown by black and green circles may be found in our recent paper [22]. The green 
circles show recent experimental data that became available to us after we completed our fit; 
Although these recent data agree reasonably well with our approximations, a refitting would 
slightly improve the agreement. 
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revised the approximations of all the integral elementary cross sections used in CEM95 [9] and 
its predecessors. We started by collecting all published experimental data from all available 
sources. Then we developed an improved, as compared with the standard Dubna INC [3], 
algorithm for approximation of cross sections and developed simple and fast approximations for 
elementary cross sections which fit very well presently available experimental data not only to 
5 GeV, the upper recommended energy for the present version of the CEM, but up to 50-100 
GeV and higher, depending on availability of data (see details in [T5l [22]). So far, we have in 
CEM03.01 new approximations for 34 different types of elementary cross sections induced by 
nucleons, pions, and gammas. Integral cross sections for other types of interactions taken into 
account in CEM03.01 are calculated from isospin considerations using the former as input. 

Examples of several compiled experimental cross sections together with our new approxi- 
mations and the old approximations from CEM95 [9j are shown in Figs. 4 and 5. We see that 
our new approximations describe indeed very well all data. Although presently we have much 
more data than 35 years ago when Barashenkov's group produced their approximations used in 
CEM95, for a number of interaction modes like the total cross sections shown in the left panel 
of Fig. 4, the original approximations also agree very well with presently available data, in the 
energy region where the Dubna INC was developed to work. This is a partial explanation of 
why the old Dubna INC [3] and the younger CEM95 [9J work so well for the majority of char- 
acteristics of nuclear reactions. On the other hand, for some modes of elementary interactions 
like the ones shown in the right panel of Fig. 4, the old approximations differ significantly from 
the present data, demonstrating the need for our recent improvements for a better description 
of all modes of nuclear reactions. 

The kinematics of two-body elementary interactions and absorption of photons and pions by 
a pair of nucleons is completely defined by a given direction of emission of one of the secondary 
particles. The cosine of the angle of emission of secondary particles in the cm. system is 
calculated by the Dubna INC [3] as a function of a random number ^, distributed uniformly in 
the interval [0,1] as 



cos e = 2^^/2 

where = M = 3, 



' N N 
n=0 n=0 



1 , (11^ 



M 

a„ = ^anfc7;'. (12) 

k=0 

The coefficients a„fc were fitted to the then available experimental data at a number of incident 
kinetic energies Tj , then interpolated and extrapolated to other energies (see details in pi^lBU] 
and references therein). The distribution of secondary particles over the azimuthal angle (p is 
assumed isotropic. For elementary interactions with more than two particles in the final state, 
the Dubna INC uses the statistical model to simulate the angles and energies of products (see 
details in [3]). 

For the improved version of the INC in CEM03.01, we use currently available experimental 
data and recently published systematics proposed by other authors and have developed new 
approximations for angular and energy distributions of particles produced in nucleon-nucleon 
and photon-proton interactions. So, for pp, np, and nn interactions at energies up to 2 GeV, 
we did not have to develop our own approximations analogous to the ones described by Eqs. 
(11) and (12), since reliable systematics have been developed recently by Cugnon et al. for 
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the Liege INC [72], then improved still further by Duarte for the BRIC code [73]; we simply 
incorporate into CEM03.01 the systematics by Duarte [73] . 

Examples of angular distributions of secondary particles from np reactions at several energies 
are shown in Fig. 6. The new approximations reproduce the experimental data much better than 
the old Dubna INC used in our previous code versions (and in several other codes developed 
from the Dubna INC). 




Figure 6: Example of twelve angular distributions of n from np elastic interactions as functions 
of '^n from 386 to 1243 MeV. The dashed lines show the old approximations from the 

Dubna INC while the solid lines are the new approximations incorporated into CEM03 and 
LAQGSM03. References to experimental data shown here by black circles may be found in our 
paper [2T] . 

In the case of 7p reactions (6) and (7), we chose another way: Instead of fitting the param- 
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eters a„ from Eq. (11) at different E-^ we found data (see, e.g., Fig. 7) and finding the energy 
dependence of parameters a^/c in Eq. (12) using the values obtained for a„, we took advantage of 
the event generator for 7p and 7n reactions from the Moscow INC ^74] kindly sent us by Dr. Igor 
Pshenichnov. That event generator includes a data file with smooth approximations through 
presently available experimental data at 50 different gamma energies from 117.65 to 6054 MeV 
(in the system where the p or n interacting with 7 is at rest) for the cm. angular distributions 
da/dQ of secondary particles as functions of G tabulated for values of G from to 180 deg., 
with the step AG = 10 deg., for 60 different channels of 7p and 'jn reactions considered by 
the Moscow INC (see details in [H]). We use part of that data file with data for reactions (6) 
and (7), and have written an algorithm to simulate unambiguously da/dQ and to choose the 
corresponding value of G for any E.y, using a single random number uniformly distributed in 
the interval [0,1]. This is straightforward due to the fact that the function ,^(cosG) 

cos 1 

^(cosG) = J da/dndcosQ^ j da/dVldcosQ 
-1 -1 

is a smooth monotonic function increasing from to 1 as cosG varies from -1 to 1. Naturally, 
when differs from the values tabulated in the data file, we perform first the needed inter- 
polation in energy. We use this procedure to describe in CEM03.01 angular distributions of 
secondary particles from reactions (6) and (7), as well as for isotopically symmetric reactions 
7 + n — i> n + 7r° and 7 + n — >■ p + vr" . 

Examples of eight angular distributions of vr"'" from 7p vr^n as functions of Qcm s 
presented in Fig. 7. We see that the approximations developed in CEM03.01 (solid histograms) 
agree much better with the available experimental data than the old Dubna INC approximations 
(11)-(12) used in all predecessors of CEM03 (dashed histograms). 

The analysis of experimental data has shown that the channel (8) of two-pion photo- 
production proceeds mainly through the decay of the isobar listed in the last Review 
of Particle Physics by the Particle Data Group as having the mass M = 1232 MeV 

7 + p — > A^^ + 7r~ , 

A++ p + n+ , (13) 

whereas the production cross section of other isobar components (|, |) are small and can be 
neglected. The Dubna INC uses the Lindenbaum-Sternheimer resonance model [75] to simulate 
the reaction (13). In this model, the mass of the isobar M is determined from the distribution 

dW 

^r.F{E,M)a{M) , (14) 

where E is the total energy of the system, F is the two-body phase space of the isobar and vr" 
meson, and a is the isobar production cross section which is assumed to be equal to the cross 
section for elastic vr+p scattering. 

The cm. emission angle of the isobar is approximated using Eqs. (11) and (12) with the 
coefficients a„fc listed in Tab. 3 of Ref. ^69j; isotropy of the decay of the isobar in its cm. system 
is assumed. 

In order to calculate the kinematics of the non-resonant part of the reaction (8) and the 
two remaining three-body channels (9) and (10), the Dubna INC uses the statistical model. 



14 



200 MeV Y + p -> 71* + n 



OAD67 
AAD60 
□ AL63 
VFR65 
OBE56 
<WA55 

— new 

— old 



40 



60 



100 



120 



140 



160 



180 



10° 



250 MeV y + p -> 71* + n 



77 



OFU77 
AFR65 

o FI72: 240 MeV 
V FI72: 260 MeV 

— new 

— old 



100 
Jdeg) 



120 



140 



160 



180 





300 MeV Y+ p -> ji"" + n 






















OFU77 






AFI72 


1 




nBE68 


1 ' 




VCL75 


1 
1 




OFR65 






— new 






— old 




. . . 1 . . . 








20 40 60 80 100 120 140 160 180 



10' 



a 10° 



400 MeV y + p -> ji"" + n 






~l_ 


i' OFU77 




] 


AFI72 






□ BE68 




1 

[ ■ 


VFR65 






. — new 




L 


: -- old 



1 ... 1 ... 1 ... 1 ... 1 .. . 


!■ 





1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 

20 40 60 80 100 120 140 160 



600 MeV y + p -> jt* + n 




OBE68 
ABE63 
□ FR65 
VDI60 

>DA73: 596 MeV 

<FI71 
— new 
-- old 



80 100 



■ — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — ■ 

790 MeV y + p -> ji* + n 











- O FU77 




" 1 




AFI71 




"I 




□ DA73 






i_ 
1 


VAL83 






1^ ■ 


- « EC67: 793 MeV 








— new 








— old 






1 . . 





80 100 



T . ' ' I ' ' ' I ' ' ' I ' ' ' I ' ' . I . ' ' I ' . . I . ' ' ' . . I . ' ' I ' ' ' I ' ' ' I ' ' ' I ' ' . I . ' ' I ' . . r 




Figure 7: Example of eight angular distributions of vr"^ from 7p — 7r~^n as functions of s 
at photon energies from 200 MeV to 1.52 GeV. The dashed lines show the old approximations 
used in the Dubna INC PRM while the solid lines are our new approximations incorporated 
into the CEM03 and LAQGSM03 codes. References to experimental data shown by symbols 
may be found in our recent paper [22] . 



The total energies of the two particles (pions) in the cm. system are determined from the 
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distribution 

___ ^ (ij _ - E^,)E^,E^JE , (15) 

and that of the third particle (nucleon, N) from conservation of energy. The actual simulation 
of such reactions is done as follows: Using a random number ^, we simulate in the beginning 
the energy of the first pion using 

where 

Then, we simulate the energy of the second pion Et,.^ according to Eq. (15) using the Monte- 
Carlo rejection method. The energy of the nucleon is calculated as Ej^ = E — E^,^ — Ej.^, 
following which we check that the "triangle law" for momenta 



IPtti - Pvra I <PN < IPtti + P 



7r2 I 



is fulfilled, otherwise this sampling is rejected and the procedure is repeated. The angles 6 and 
(p of the pions are sampled assuming an isotropic distribution of particles in the cm. system, 

cos e^, = 2^1 - 1, COS = 2^2 - 1, V'tti = 27r^3, (f^^ = 27r^4, 

and the angles of the nucleon are defined from momentum conservation, p^r = —{Pm + P7r2)- 
More details on our new approximations for differential elementary cross sections may be found 

in [SUES]. 

The Pauli exclusion principle at the cascade stage of the reaction is handled by assuming 
that nucleons of the target occupy all the energy levels up to the Fermi energy. Each simulated 
elastic or inelastic interaction of the projectile (or of a cascade particle) with a nucleon of the 
target is considered forbidden if the "secondary" nucleons have energies smaller than the Fermi 
energy. If they do, the trajectory of the particle is traced further from the forbidden point and 
a new interaction point, a new partner and a new interaction mode are simulated for the traced 
particle, etc., until the Pauli principle is satisfied or the particle leaves the nucleus. 

In this version of the INC, the kinetic energy of the cascade particles is increased or decreased 
as they move from one of the seven potential regions (zones) to another, but their directions 
remain unchanged. That is, in our calculations, refraction or reflection of cascade nucleons 
at potential boundaries is neglected. CEM03.01 allows us to take into account refractions and 
reflections of cascade nucleons at potential boundaries; for this, one needs to change the value of 
the parameter irefrac from to 1 in the subroutine initial. But this option provides somewhat 
worse overall agreement of calculations with some experimental data, therefore the option of 
no refractions/reflections was chosen as the default in CEM03.01. 

The INC in CEM does not take into account the so-called "trawling" effect |3j. That is, in 
the beginning of the simulation of each event, the nuclear density distributions for the protons 
and neutrons of the target are calculated according to Eq. (1) and a subsequent decrease of the 
nuclear density with the emission of cascade particles is not taken into account. Our detailed 
analysis of different characteristics of nucleon- and pion-induced reactions for targets from C 
to Am has shown that this effect may be neglected at incident energies below about 5 GeV in 
the case of heavy targets like actinides and below about 1 GeV for light targets like carbon. 
At higher incident energies the progressive decrease of nuclear density with the development of 



16 



the intranuclear cascade has a strong influence on the calculated characteristics and this effect 
has to be taken into account [3]. Therefore, in transport codes that use as event generators 
both CEM03.01 [23] and our high-energy code LAQGSM03.01 [24], we recommend simulating 
nuclear reactions with CEM03.01 at incident energies up to about 1 GeV for light nuclei like C 
and up to about 5 GeV for actinide nuclei, and to switch to simulations using LAQGSM03.01, 
which considers the "trawling" effect, at higher energies of transported particles. 

An important ingredient of the GEM is the criterion for transition from the intranuclear 
cascade to the preequilibrium model. In conventional cascade-evaporation models (like the 
Bertini ING used in MGNPX fast particles are traced down to some minimal energy, the 
cutoff energy Tcut (or one compares the duration of the cascade stage of a reaction with a cutoff 
time, in "time-like" ING models, such as the Liege ING [T2j). This cutoff is usually less than 
~ 10 MeV above the Fermi energy, below which particles are considered to be absorbed by the 
nucleus. The GEM uses a different criterion to decide when a primary particle is considered to 
have left the cascade. 

An effective local optical absorptive potential Wopt.mod.i''^) is defined from the local interac- 
tion cross section of the particle, including Pauli-blocking effects. This imaginary potential is 
compared to one defined by a phenomenological global optical model Wopt. exp.i''^) ■ charac- 
terize the degree of similarity or difference of these imaginary potentials by the parameter 

V =1 (Wopt.mod. — Wopt. exp.) /Wopt. exp. \ ■ (16) 

When V increases above an empirically chosen value, the particle leaves the cascade, and is 
then considered to be an exciton. From a physical point of view, such a smooth transition from 
the cascade stage of the reaction seems to be more attractive than the "sharp cut-off" method. 
In addition, as was shown in Ref. |2j, this improves the agreement between the calculated and 
experimental spectra of secondary nucleons, especially at low incident energies and backward 
angles of the detected nucleons (see e.g., Figs. 3 and 11 of Ref. P]). More details about this 
can be found in [21 [TTl [76] . 

GEM03.01 uses a fixed value V = 0.3 (at incident energies below 100 MeV), just as all 
its predecessors did. With this value, we find that the cascade stage of the GEM is generally 
shorter than that in other cascade models. This fact leads to an overestimation of preequi- 
librium particle emission at incident energies above about 150 MeV, and correspondingly to 
an underestimation of neutron production from such reactions, as was established in Ref. |17j . 
In Ref. [T7j, this problem was solved temporarily in a very rough way by using the transition 
from the ING to the preequilibrium stage according to Eq. (16) when the incident energy of 
the projectile is below 150 MeV, and by using the "sharp cut-off" method with a cutoff energy 
Tcut = 1 MeV for higher incident energies. This "ad hoc" rough criterion solved the problem of 
underestimating neutron production at high energies, providing meanwhile a reasonably good 
description of reactions below 150 MeV. But it provides an unphysical discontinuity in some 
observables calculated by MGNPX using GEM2k [T7] as an event generator, observed but not 
understood by Broeders and Konobeev [77]. In GEM03.01, this problem is solved by using a 
smooth transition from the first criterion to the second one in the energy interval from 75 to 
225 MeV, so that no discontinuities are produced in results from GEM03.01. 

Beside the changes to the Dubna ING mentioned above, we also made in the ING a number of 
other improvements and refinements, such as imposing momentum-energy conservation for each 
simulated event (the Monte-Garlo algorithm previously used in the GEM provided momentum- 
energy conservation only statistically, on the average, but not exactly for each simulated event) 
and using real binding energies for nucleons in the cascade instead of the approximation of 
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a constant separation energy of 7 MeV used in previous versions of the CEM. We have also 
improved many algorithms used in the Monte-Carlo simulations in many subroutines, decreasing 
the computing time by up to a factor of 6 for heavy targets, which is very important when 
performing practical simulations with transport codes like MCNPX or MARS. 

Let us mention that in the CEM the initial configuration for the preequilibrium decay 
(number of excited particles and holes, i.e. excitons no = po + ^o, excitation energy Eq, lin- 
ear momentum Pq, and angular momentum Lq of the nucleus) differs significantly from that 
usually postulated in exciton models. Our calculations [21 [THl [79] have shown that the dis- 
tributions of residual nuclei remaining after the cascade stage of the reaction, i.e. before the 
preequihbrium emission, with respect to no, Po; ^O; -E'o; -F*0; and Lq are rather broad. ^ 

CEM03.01 (just like LAQGSM03.01 and many other INC-based models) calculates the total 
reaction cross section, cxin, by the Monte-Carlo method using the geometrical cross section, 
Cgeom, and the number of inelastic, A^j„, and elastic, N^u simulated events, namely: = 
(^geomNin/ {Nin + N^i). The valuc of the total reaction cross section calculated this way is printed 
in the beginning of the CEM03.01 output labeled as Monte Carlo inelastic cross section. This 
approach provides a good agreement with available data for reactions induced by nucleons, 
pions, and photons at incident energies above about 100 MeV, but is not reliable enough at 
energies below 100 MeV (see, e.g., Figs. 8 and 9 below). 

To address this problem, we have incorporated [11] into CEM03.01 the NASA systematics by 
Tripathi et al. [HO] for all incident protons and for neutrons with energies above the maximum in 
the NASA reaction cross sections, and the Kalbach systematics [81] for neutrons of lower energy. 
For reactions induced by monochromatic and bremsstrahlung photons, we incorporate [22] into 
CEM03.01 the recent systematics by Kossov [82]. Details on these systematics together with 
examples of several total inelastic cross sections calculated with them compared with available 
experimental data may be found in [T^l [22] • Our analysis of many different reactions has 
shown that at incident energies below about 100 MeV these systematics generally describe the 
total inelastic cross sections better that the Monte-Carlo method does, and no worse than the 
Monte-Carlo method at higher energies. Therefore we choose these systematics as the default 
for normalization of all CEM03.01 results. The total reaction cross sections calculated by these 
systematics are printed in the CEM03.01 output labeled as Inelastic cross section used here. 
(Of course, users may re-normalize all the CEM03.01 results to the Monte-Carlo total reaction 
cross sections by making a small change to the code in the subroutine typeout). 

Let us note, however, that in applications, when CEM03.01 and LAQGSM03.01, or any other 
codes, are used as event generators in transport codes, it does not matter how they calculate 
the total reaction cross sections (normalization): All transport codes use their own routines or 
systematics to calculate the total elastic and inelastic cross sections of the projectiles, before 
starting to simulate with an event generator an inelastic interaction of the traced projectile 
with a nucleus of the thick target. 

To summarize this Section, in comparison with the initial version of the Dubna INC [3l [1] 

^Unfortunately, this fact was misunderstood by the authors of the code HETC-3STEP [59 . In spite of 
the fact that it has been stressed exphcitly, and figures with distributions of excited nuclei after the cascade 
stage of a reaction with respect to the number of excitons and other characteristics were shown in a number 
of publications (see, e.g., Fig. 5 in Ref. [2], Fig. 1 in Ref. [79 , p. 109 in Ref. [78], and p. 706 in Ref. 26J), the 
authors of Ref. |59| misstated this fact as "Gudima et al. assumed the state of two particles and one hole at 
the beginning ■ ■ ■ Hence, their assumption is not valid for the wide range of incident energy", claiming this as 
a weakness of the CEM and a priority of the code HETC-3STEP, where smooth distributions of excited nuclei 
after the cascade stage of reactions with respect to ng are used. This had already been done in the CEM [1] [5] . 
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Figure 8: Total reaction cross sections for n- and p-induced reactions on Al calculated by 
CEM03.01 and its predecessor CEM97 with experimental data compiled by Barashenkov [83] 
and calculations from the HMS-ALICE code [Ml. 
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Figure 9: Examples of total photoabsorption cross sections for ^^'^Sn, ^^^Ta, ^^''Au, ^^'^Pb, 
^^^Th, and ^ss-y functions of photon energy. The red lines marked as "GABS. FOR" are 
results by a CEM03.01 subroutine written to reproduce Kossov's [HS] systematics, as described 
in [22]. The green line marked as "LANL", "KAERI", or "BOFODM" show the evaluations 
by LANL, KAERI, or by a collaboration between IPPE/Obninsk and CDFE/Moscow (the 
BOFOD(MOD) Library) from the IAEA Photonuclear Data Library [85]. Results from the 
photonuclear version of CEM95 [7\ and from CEM2k as modified for MCNPX by Gallmeier 
[86] are shown by the blue and brown dashed lines, respectively. References to experimental 
data shown by different symbols may be found in our recent paper [22] . 
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used in CEM95 [9], for CEM03.03 we have: 

1) developed better approximations for the total elementary cross sections; 

2) developed new approximations to describe more accurately experimental elementary en- 
ergy and angular distributions of secondary particles from hadron-hadron and photon-hadron 
interactions; 

3) normalized the photonuclear reactions to detailed systematics developed by M. Kossov 
and the nucleon-induced reactions, to NASA and Kalbach systematics; 

4) changed the condition for transition from the INC stage of a reaction to preequilibrium; 
on the whole, the INC stage in CEM03.01 is longer while the preequilibrium stage is shorter in 
comparison with previous versions; 

5) incorporated real binding energies for nucleons in the cascade instead of the approxima- 
tion of a constant separation energy of 7 MeV used in the initial versions of the CEM; 

6) imposed momentum-energy conservation for each simulated even (provided only "on the 
average" by the initial versions); 

7) changed and improved the algorithms of many INC routines and almost all INC routines 
were rewritten, which speeded up the code significantly; 

8) fixed some preexisting bugs. 

On the whole, the INC of CEM03.01 describes nuclear reactions better and much faster 
than the the initial version of the Dubna INC [3],|1] used in CEM95 [9J. One example of results 
by the INC from our CEM03.01 is shown in Fig. 10, namely, 7r° spectra from 500 MeV vr^ + Cu. 




100 200 300 400 500 

T (MeV) 

Figure 10: Experimental 7r° spectra from 500 MeV vr" + Cu [271 EZ] compared with CEM03.01 
results. Let us recall here that as pions are produced by CEM03.01 only at the INC stage of 
reactions, calculated pion spectra do not depend at all on how other reaction mechanisms like 
coalescence, evaporation, fission, or Fermi breakup are calculated. 



21 



3.2. The INC of LAQGSM03.03 



The INC of LAQGSM03.03 is described with a recently improved version [2H [88] of 
the time-dependent intranuclear cascade model developed initially at JINR in Dubna, often 
referred to in the literature as the Dubna intranuclear Cascade Model, DCM (see [BH] and 
references therein). The DCM models interactions of fast cascade particles ("participants") 
with nucleon spectators of both the target and projectile nuclei and includes as well interactions 
of two participants (cascade particles). It uses experimental cross sections at energies below 4.5 
GeV/nucleon, or those calculated by the Quark-Gluon String Model ^0j-[96j at higher energies 
to simulate angular and energy distributions of cascade particles, and also considers the Pauli 
exclusion principle. 

In contrast to the CEM03.01 version of the INC described above, DCM uses a continuous 
nuclear density distribution (instead of the approximation of several concentric zones, where 
inside each the nuclear density is considered to be constant); therefore, it does not need to 
consider refraction and reflection of cascade particles inside or on the border of a nucleus. 
It also keeps track of the time of an intranuclear collision and of the depletion of the nuclear 
density during the development of the cascade (the so-called "trawling effect" mentioned above) 
and takes into account the hadron formation time (see Fig. 11). 




Before starting to simulate an INC event, position of The projectile interacts (in point A) with the nearest 
all IntraNudear nucleons are simulated and "frozen" target nucleon met Inside the cylinder with the radius r 




*'i(2.3,,,,) is the formation time of the cascade particle #1 (2,3,...) 
If tj < t,, t2 < tj,..., and t2 > t'j, particle #2 interacts first in point C 

Figure 11: An illustrative scheme of a target nucleus, of interaction points of cascade particles 
(participants) with intranuclear nucleons (spectators), and of selection of the corresponding 
time of such interactions, as performed in the INC used in LAQGSM. 

In the INC used in LAQGSM, all the new approximations developed recently for the INC 
of CEM03.01 to describe total cross sections and elementary energy and angular distributions 
of secondary particles from hadron-hadron interactions have been incorporated [21]. Then, a 
new high-energy photonuclear reaction model based on the of the event generators for 7p and 
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reactions from the Moscow INC [71] kindly provided us by Dr. Igor Pshenichnov, and on 
the latest photonuclear version of CEM |22] was developed and incorporated [88] into the INC 
of LAQGSM, which allows us to calculate reactions induced by photons with energies of up to 
tens of GeV. Finally, the algorithms of many LAQGSM INC routines were changed and some 
INC routines were rewritten, which speeded up the code significantly; some preexisting bugs in 
the DCM were fixed; many useful comments were added [T2] . 

In the latest version of LAQGSM, LAQGSM03.03 [I2], the INC was modified for a better 
description of nuclear reactions at very high energies (above 20 GeV/nucleon), namely: 

1) The latest fits to currently available evaluated experimental databases for the total and 
elastic tt~^p, vr~p, pp, and pn cross sections (see Chapter 40 in the last Review of Particle Physics 
[98] and references therein) have been incorporated into LAQGSM. LAQGSM03.03 uses now 
these approximations at energies above 20-30 GeV, and its own approximations developed for 
CEM03.01 [23J at lower energies. 

2) Previously, LAQGSM was used only at energies below 800 GeV. In [12], the possibility 
of using LAQGSM03.03 at ultra-relativistic energies, above 1 TeV was studied. It was found 
that to describe ultra-high energy reactions, the value of the parameter a± = 0.51 GeV/c in 
the transverse momentum distribution of the constituent quarks of QGSM (see Eq. (12) in ^3] 
or Eq. (10) in Ref. [93]) has to be increased from 0.51 GeV/c at Tp < 200 GeV [ISj to ^ 2 
GeV/c at Tp ~ 21 TeV. 

More details on the INC and other nuclear reaction models considered by different versions 
of LAQGSM may be found in Refs. [121 1131 [2H EHl EZ] • Several examples of recent results by 
LAQGSM are shown in Figs. 12-16. 

Fig. 12 shows a test of LAQGSM03.03 on inclusive pion production spectra in proton- 
beryllium collisions at 6.4, 12.3, and 17.5 GeV/c obtained from data taken by the already quite 
old E910 measurements at Brookhaven National Laboratory, but analyzed and published only 
several months ago [99] . Let us recall again that pions are produced only at the INC stage of 
reactions, without any contributions from other reaction mechanisms, so that such results test 
just the INC part of any model. We see that LAQGSM03.03 describes these pion spectra quite 
well, just as we obtained and published with previous versions of LAQGSM for other spectra 
of different particles measured by the E910 experiment. 

Fig. 13 presents part of the recent extensive experimental data on fragmentation cross 
sections of ^®Si on H, C, Al, Cu, Sn, and Pb at energies from 290 to 1200 MeV/nucleon [TOO] . 
Such measurements are of interest for NASA to plan long-duration space-flights and to test 
the models used to evaluate radiation exposure in flight, and were performed at many incident 
energies in this energy range at the Heavy Ion Medical Accelerator in Chiba (HIMAC) and at 
Brookhaven National Laboratory (see details in jlOO] and references therein). We calculate in 
our model practically all these data, but here limit ourselves to examples of results for only 
three energies, for each measured target. For comparison, we present in Fig. 13 results from 
both LAQGSM03.03 (solid lines) and its predecessor LAQGSM03.02 (dashed lines). In general, 
LAQGSM03.03 describes these new data slightly better than LAQGSM03.02 [63], although this 
is not obvious on the scale of the figure. The agreement of our calculations with these data is 
excellent, especially considering that the results presented in this figure, just as all our other 
results shown in these lectures, are obtained without fitting any parameters in our codes; 
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Figure 12: Measured inclusive forward 7r+ and vr spectra from 6.4, 12.3, and 17.5 GeV/c p + 
^Be |99j compared with LAQGSM03.03 results at angles of detection as indicated in the plots. 



For reactions induced by 6.4 GeV/c protons, we also show LAQGSM03.03 predictions 
unmeasured spectra at 90 and 159 degrees. 
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Figure 13: Atomic-number dependence of the fragment-production cross sections from the 
interactions of ^^Si of about 270, 560, and 1150 MeV/nucleon with H, C, Al, Cu, Sn, and Pb, 
as indicated. Filled circles are measurements by Zeitlin et al. [100]; solid lines are results from 
LAQGSM03.03 [I2], while dashed lines are results from LAQGSM03.02 
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we simply input A and Z of the projectile and target and the energy of the projectile, then 
calculate without changing or fitting anything. 

Fig. 14. shows proton spectra spectra at 30, 60, 90, 120 and 150 degrees from interaction 
of bremsstrahlung 7 quanta of maximum energy 4.5 GeV with ^^C, ^''Al, ^^Cu, and ^°^Pb. 
Experimental data shown by symbols in the figures are quite old, measured about 30 years 
ago by Alanakyan et al. |101l 1102] , however, to the best of our knowledge, we have described 
with LAQGSM03.01 these data for the first time: We do not know of any publications or 
oral presentations where these measurements were reproduced by a theoretical model, event 
generator, or transport code. 
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Figure 14: Proton spectra at 30, 60, 90, 120 and 150 degrees from interaction of bremsstrahlung 
7 quanta of maximum energy 4.5 GeV with ^^C, ^^Al, ^^Cu, and ^°^Pb. Experimental values 
shown by symbols are from |101l I102j while histograms show results from LAQGSM03.01. 
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Fig. 15. shows an example of neutron spectra from intermediate-energy nucleus-nucleus 
interactions studied lately by many authors because of a great interest in such reactions for 
radiotherapy and because such processes contribute substantially to the dose and dose equiv- 
alent in space-flight. Namely, in Fig. 15, we compare results from LAQGSM03.01 for neutron 
spectra at 5, 10, 20, 30, 40, 60, and 80 degrees from interactions of 600 MeV/nucleon ^°Ne, on 
^^C, ^^Al, ^"^Cu, and ^''^Pb with measurements by Iwata et al. |103j and calculations by JQMD 
1 104] and HIC [105]. We see that LAQGSM describes these data reasonably well, generally as 
well as or better than do JQMD or HIC. 



\ h I 1 1 1 1 




10 10 10 10 10 10 

Tn (MeV) 

Figure 15: Comparison of measured |103j double differential cross sections of neutrons from 
600 MeV/nucleon ^ONe, on ^^c, 27ai, 64q^^ 208pb ^j^j^ LAQGSM03.01 resuhs and 
calculations by JQMD |104] and HIC |105j . Experimental data for these reactions on Al are 
not yet available so we present here |106j only our predictions from LAQGSM03.01. 

Finally, Fig. 16 shows an example of product mass yields measured recently at GSI in 
inverse kinematics |107[ llOSj . namely, product yields of 9 isotopes from Zn to Hg produced 
from interactions of a ^^^U beam with a liquid-hydrogen target as calculated by LAQGSM03.01 
as a stand-alone code and by the transport code MARS15 pT] using LAQGSM03.01 as its 
event generator |109] . The results from MARS15 using LAQGSM03.01 agree very well with 
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the results from LAQGSM03.01 as a stand-alone code and with experimental data. 



1 GeV p + ''°U 
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Figure 16: Experimental mass yield in 1 GeV p + ^ss^j measured in inverse kinematics at 
GSI |107l 1108] compared with results by MARS15 using LAQGSM03.01 as its event generator 
(dashed lines) and by LAQGSM03.01 as a stand-alone code (solid lines). 



4. The Coalescence Model 

When the cascade stage of a reaction is completed, CEMOS.Ox and LAQGSMOS.Ox use 
the coalescence model described in Refs. |89l IllOj to "create" high-energy d, t, ^He, and ''He 
by final-state interactions among emitted cascade nucleons, already outside of the target nu- 
cleus. In contrast to most other coalescence models for heavy-ion-induced reactions, where 
complex-particle spectra are estimated simply by convolving the measured or calculated inclu- 
sive spectra of nucleons with corresponding fitted coefficients (see, e.g., [lllj and references 
therein), CEMOS.Ox and LAQGSMOS.Ox use in their simulations of particle coalescence real 
information about all emitted cascade nucleons and do not use integrated spectra. We assume 
that all the cascade nucleons having differences in their momenta smaller than pc and the cor- 
rect isotopic content form an appropriate composite particle. This means that the formation 
probability for, e.g. a deuteron is 

Wd{p, = j j dppdpnP^iPp, h)p^{pn, b)5{pp + Pn- p)Q{Pc - \Vp - Pn\) , (17) 

where the particle density in momentum space is related to the one-particle distribution function 
/by 

p^(p,6) = I dff{f,p,b). (18) 
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Here, b is the impact parameter for the projectile interacting with the target nucleus and the 
superscript index C shows that only cascade nucleons are taken into account for the coalescence 
process. The coalescence radii Pc were fitted for each composite particle in Ref. [89j to describe 
available data for the reaction Ne+U at 1.04 GeV/nucleon, but the fitted values turned out to 
be quite universal and were subsequently found to satisfactorily describe high-energy complex- 
particle production for a variety of reactions induced both by particles and nuclei at incident 
energies up to about 200 GeV/nucleon, when describing nuclear reactions with different versions 
of LAQGSM [121 [631 [13] or with its predecessor, the Quark-Gluon String Model (QGSM) [901 • 
These parameters are: 

Pc{d) = 90 MeV/c; Pc{t) = PcC^e) = 108 MeV/c; PcCHe) = 115 MeV/c . (19) 

As the INC of CEM03.0x is different from those of LAQGSM or QGSM, it is natural to expect 
different best values for p^ as well. Our recent studies show that the values of parameters Pc 
defined by Eq. (19) are also good for CEM03.01 for projectile particles with kinetic energies Tq 
lower than 300 MeV and equal to or above 1 GeV. For incident energies in the interval 300 MeV 
< Tq < 1 GeV, a better overall agreement with the available experimental data is obtained 
by using values of equal to 150, 175, and 175 MeV/c for d, t(^He), and ^He, respectively. 
These values of Pc are fixed as defaults in CEM03.01. If several cascade nucleons are chosen to 
coalesce into composite particles, they are removed from the distributions of nucleons and do 
not contribute further to such nucleon characteristics as spectra, multiplicities, etc. 

In comparison with the initial version [89l[lT0], in CEM03.0x and LAQGSM03.0x, several 
coalescence routines have been changed/deleted and have been tested against a large variety of 
measured data on nucleon- and nucleus-induced reactions at different incident energies. 

Two examples of results from the coalescence model are shown in Figs. 17 and 18. We see 
that for a reaction between an intermediate-energy neutron and a medium-mass nucleus (Fig. 
17; 96 MeV n + Fe), where the mean multiplicity of the secondary nucleons is small, the con- 
tribution from coalescence to the total angle-integrated energy spectra of complex particles is 
very low, less than a few percent. On the other hand, in the case of interaction of a high-energy 
proton with a heavy nucleus (Fig. 18, 70 GeV p + Pb), very energetic secondary complex 
particles produced at forward angles: deuterons with 14 GeV/c and tritons with 19 GeV/c 
were measured |113j at 160 mrad. Probably all such extremely energetic d and t at forward 
angles are produced only via coalescence of complex particles particles from energetic nucleons 
emitted during the INC stage of the reaction: We do not know any other interaction mecha- 
nisms that would produce d and t of such high energies from this reaction. The coalescence 
model in LAQGSM03.0x reproduces these experimental spectra quite well. It is clear that the 
coalescence mechanism is more important for high-energy heavy-ion reactions, where the mul- 
tiplicity of secondary INC nucleons is much higher than in the case of nucleon-induced reactions. 

5. Preequilibrium Reactions 

The subsequent preequilibrium interaction stage of nuclear reactions is considered by our 
current GEM and LAQGSM in the framework of the latest version of the Modified Exciton 
Model (MEM) [3 |5] as implemented in CEM03.01 [23]. At the preequilibrium stage of a 
reaction we take into account all possible nuclear transitions changing the number of excitons 
n with An = +2, —2, and 0, as well as all possible multiple subsequent emissions of n, p, d, 
t, ^He, and ^He. The corresponding system of master equations describing the behavior of a 
nucleus at the preequilibrium stage is solved by the Monte-Carlo technique [Il[2]. 
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Figure 17: Angle-integrated energy spectra of n, p, d, t, ^B.e, and ^He from 96 MeV n + 
^^Fe calculated by CEM03.01 compared with experimental data by Blideanu et al. |112j . The 
black histograms show the total calculated spectra (all reaction mechanisms considered by 
CEMOS.Ox), while the color histograms show separately contributions to the total spectra from 
INC, coalescence, preequilibrium, and evaporation, respectively, as indicated in the correspond- 
ing legends of plots. 



30 




Figure 18: Invariant cross sections EcPa/d^p for forward production of p, d, and t at 160 mrad 
(9.17 deg) as functions of particle momentum p from 70 GeV protons on ^°^Pb. Experimental 
data shown by stars are from Ref. [113] , while calculations by LAQGSM03.01, LAHET3 ^114j . 
and MARS 14 [51] are shown by histograms, as indicated in the legends. 



For a preequilibrium nucleus with excitation energy E and number of excitons n = p + h, 
the partial transition probabilities changing the exciton number by An are 

\An{p,h,E) = —\MAn?UAn{p,KE) . (20) 

n 

The emission rate of a nucleon of type j into the continuum is estimated according to the 
detailed balance principle 



T,{p,h,E) = j Xi{p,h,E,T)dT 
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KiP,KE,T) = ^-^^^,^,ip,h f^P 'u^Ie)' ^^ ^^-(^)' (21) 

where Sj, Bj, V^, and /Zj are the spin, binding energy, Coulomb barrier, and reduced mass 
of the emitted particle, respectively. The factor 3?j(j9, h) ensures the condition for the exciton 
chosen to be the particle of type j and can easily be calculated by the Monte-Carlo technique. 

Assuming an equidistant level scheme with the single-particle density g, we have the level 
density of the n-exciton state as [115] 

coip, h, E) = ,f,7 \ . (22) 

This expression should be substituted into Eq. (21). For the transition rates (20), one needs 
the number of states taking into account the selection rules for intranuclear exciton-exciton 
scattering. The appropriate formulae have been derived by Williams |116] and later corrected 
for the exclusion principle and indistinguishability of identical excitons in Refs. |117l 1118] : 

1 [gE-A{p+l,h + l)f^gE-A{p+l,h+iy-^ 
uj+{p,h,E) = -g- 



T n + 1 L gE-A{p,h) 

u;o{p,h,E) = ^g\^^—^iPM[p(p-l)+4ph + h{h-l)], 

uj_{p,h,E) = ^gph{n-2) , (23) 

where A{p, h) = {p'^ + h'^ +p — h)/A — h/2. By neglecting the difference of matrix elements with 
different An, M+ = M_ = Mq = M, we estimate the value of M for a given nuclear state by 
associating the A+(p, h, E) transition with the probability for quasi- free scattering of a nucleon 
above the Fermi level on a nucleon of the target nucleus. Therefore, we have 

< (T{Vrel)Vrel > _ TT ^^^^ (^[(^E - ^(p + 1, /l + 1)] V gE - Ajp + I, h + l) l "'^ 

Vint ^' ' n + 1 I gE-Aip,h) J ■ ^ ' 

Here, Vmt is the interaction volume estimated as Vint = |7r(2rc + A/27r)^, with the de Broglie 
wave length A/27r corresponding to the relative velocity Vrei = ^/^Trei/rriN- A value of the 
order of the nucleon radius is used for Tc in the CEM: rc = 0.6 fm. 

The averaging in the left-hand side of Eq. (24) is carried out over all excited states taking 
into account the Pauli principle in the approximation 

< a{Vrel)Vrel >-< (^{Vrel) >< Vrel > ■ (25) 

The averaged cross section < a{vrei) > is calculated by the Monte-Carlo simulation method 
and by introducing a factor rj effectively taking into account the Pauli principle exactly as is 
done in the Fermi-gas model (see, e.g., |119j )^ 

(7{Vrel) = ^Wpp{Vrel) + (Tpn{Vrel)]v{TF /T) , whcrC (26) 

{ \ — f ^ ~ f-^' if x < 0.5 , , . 

^W = i l_Za;+ 2^(2-1)5/2 if a; > 0.5. ^^^j 



•^Unfortunately, formula (27) as presented in Ref. [2] had some misprints; in the prior publication [T], it was 
correct. 
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Here, Vrei is the relative velocity of the excited nucleon (exciton) and the target nucleon in units 
of the speed of light and T is the kinetic energy of the exciton. The free-particle interaction 
cross sections (Jpp{vrei) and apn{vrei) in Eq- (26) are estimated using the relations suggested by 
Metropolis et al. [li] 

, , 10.63 29.92 

0-pp{Vrel) = h 42.9 , 

Kel ^rel 

, , 34.10 82.2 , , 

OpniVre^ = — + 82.2 , 28 

Kel Vrel 

where the cross sections are given in mb. 

The relative kinetic energy of colliding particles necessary to calculate < Vrei > and the 
factor rj in Eqs. (26,27) are estimated in the so-called "right-angle-collision" approximation [5], 
i.e. as a sum of the mean kinetic energy of an excited particle (exciton) measured from the 
bottom of the potential well Tp = Tp + E/n plus the mean kinetic energy of an intranuclear 
nucleon partner T/v = 3Tp/5, that is Trei = Tp + T/v = 8Tp/5 + E/n. 

Combining (20), (22) and (24), we get finally for the transition rates: 



X+{P,h,E) 
Ao(p, h,E) 



< a{Vrel)Vrel > 



int 



< (T{Vrel)Vrel > ^ + 1 

Vint n 



gE-A{p,h) ]"+ip(p- 1) +4p/i + /i(/i - 1) 



gE-A{p + l,h + l)i gE-A{p,h) 



. ( < (^{vrei)vrei > \ gE-A{p,h) ] "+i p/i(n + 1) (n - 2) 

""-^P^^^^^ - - [gE-A{p+l,h + l)\ [gE-Aip,hW ' ^ ^ 



int 



CEM considers the possibility of fast d, t, ^He, and ^He emission at the preequilibrium 
stage of a reaction in addition to the emission of nucleons. We assume that in the course 
of a reaction pj excited nucleons (excitons) are able to condense with probability jj forming 
a complex particle which can be emitted during the preequilibrium state. A modification of 
Eq. (21) for the complex-particle emission rates is described in detail in Refs. [U [2]. The 
"condensation" probability is estimated in those references as the overlap integral of the 
wave function of independent nucleons with that of the complex particle (cluster) 

7. ^ P%v,/vy^-' = p%p,/Ay^-' . (30) 

This is a rather crude estimate. In the usual way the values jj are taken from fitting the 
theoretical preequilibrium spectra to the experimental ones, which gives rise to an additional, as 
compared to (30), dependence of the factor •yj onpj and excitation energy (see, e.g., Refs. |12H 
I122j ). for each considered reaction. 

The single-particle density gj for complex-particle states is found in the CEM by assuming 
the complex particles move freely in a uniform potential well whose depth is equal to the binding 
energy of this particle in a nucleus [2] 

+ (31) 

As we stated previously, this is a crude approximation and it does not provide a good 
prediction of emission of preequilibrium a particles (see, e.g., [TB] and references therein). In 
CEM03.0x, to improve the description of preequilibrium complex-particle emission, we estimate 
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7j by multiplying the estimate provided by Eq. (30) by an empirical coefficient Mj{A, Z,To) 
whose values are fitted to available nucleon-induced experimental complex-particle spectra. 
We fix the fitted values of Mj{A, Z,Tq) in data commons of CEM03.0x and complement them 
with routines gambetn and gambetp for their interpolation outside the region covered by 
our fitting. As shown already above in Fig. 17 and proved in several more figures below, 
after fitting Mj{A, Z,To), CEM and LAQGSM describe quite well the measured spectra of all 
complex particles, providing a much better agreement with experimental data than all their 
predecessors did. 

CEM and LAQGSM predict forward-peaked (in the laboratory system) angular distributions 
for preequilibrium particles. For instance, CEM03.0x assumes that a nuclear state with a given 
excitation energy E* should be specified not only by the exciton number n but also by the 
momentum direction Q. Following Ref. [123] . the master equation (11) from Ref. [2] can be 
generalized for this case provided that the angular dependence for the transition rates A+, Aq, 
and A_ (Eq. (29) is factorized. In accordance with Eqs. (24) and (25), in the CEM it is assumed 
that 

< a >^< a > F{n) , (32) 

where 

^^^^ - JdQ'dafr^^/dQ' ■ ^^^^ 

The scattering cross section da^^'^'^/dQ is assumed to be isotropic in the reference frame of the 
interacting excitons, thus resulting in an asymmetry in both the nucleus center-of-mass and 
laboratory frames. The angular distributions of preequilibrium complex particles are assumed 
[2] to be similar to those for the nucleons in each nuclear state. 

This calculation scheme is easily realized by the Monte-Carlo technique. It provides a good 
description of double differential spectra of preequilibrium nucleons and a not-so-good but still 
satisfactory description of complex-particle spectra from different types of nuclear reactions 
at incident energies from tens of MeV to several GeV. For incident energies below about 200 
MeV, Kalbach [124] has developed a phenomenological systematics for preequilibrium-particle 
angular distributions by fitting available measured spectra of nucleons and complex particles. 
As the Kalbach systematics are based on measured spectra, they describe very well the double- 
differential spectra of preequilibrium particles and generally provide a better agreement of 
calculated preequilibrium complex-particle spectra with data than does the CEM approach 
based on Eqs. (32,33). This is why we have incorporated into CEM03.0x and LAQGSM03.0x 
the Kalbach systematics [124] to describe angular distributions of both preequilibrium nucleons 
and complex particles at incident energies up to 210 MeV. At higher energies, we use the CEM 
approach based on Eqs. (32,33). 

By "preequilibrium particles" we mean particles which are emitted after the cascade stage 
of a reaction but before achieving statistical equilibrium at a time teq, which is fixed by the 
condition A+(neq, E) = \_{n^q, E) from which we get 



"'eq — 



. (34) 



At t > teq (or n > rieq), the behavior of the remaining excited compound nucleus is described 
in the framework of both the Weisskopf-Ewing statistical theory of particle evaporation [125] 
and fission competition according to Bohr- Wheeler theory [126J. 

The parameter g entering into Eqs. (29) and (34) is related to the level-density parameter 
of single-particle states a = tt'^q/G. At the preequilibrium stage, we calculate the level-density 
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parameter a with our own approximation [16] in the form proposed initially by Ignatyuk et al. 
|127j . following the method by lljinov et al. |128] : 



a(Z, N,E*) 



~a{A) 



l + 5Wgs{Z,N) 



f(E* - A) 



(35) 



E* - A 



where 



d{A) 



aA + (3A^/^B, 



(36) 



is the asymptotic Fermi-gas value of the level-density parameter at high excitation energies. 
Here, Bg is the ratio of the surface area of the nucleus to the surface area of a sphere of the 
same volume (for the ground state of a nucleus, ~ 1), and 



E* is the total excitation energy of the nucleus, related to the "thermal" energy U by: U = 
E* — En — A, where Er and A are the rotational and pairing energies, respectively. 

We use the shell correction 6Wgs{Z, N) by MoUer et al. |129j and the pairing energy shifts 
from MoUer, Nix, and Kratz [130] . The values of the parameters a, j3, and 7 were derived in 
Ref. [16] by fitting the the same data analyzed by lljinov et al. [128] (we discovered that lljinov 
et al. used 11 /^/A for the pairing energies A in deriving their level-density systematics instead 
of the value of 11l\l~A stated in Ref. |128] and we also found several misprints in the nuclear 
level-density data shown in their Tables. 1 and 2 used in the fit). We find: 



As mentioned in Section 3.1, the standard version of the CEM ^ provides an overestimation 
of preequilibrium particle emission from different reactions we have analyzed (see more details 
in [T71 [IB])- One way to solve this problem suggested in Ref. [TT] is to change the criterion for 
the transition from the cascade stage to the preequilibrium one, as described in Section 3.1. 
Another easy way suggested in Ref. [T7] to shorten the preequilibrium stage of a reaction is 
to arbitrarily allow only transitions that increase the number of excitons. An = +2, z.e., only 
allow the evolution of a nucleus toward the compound nucleus. In this case, the time of the 
equilibration will be shorter and fewer preequilibrium particles will be emitted, leaving more 
excitation energy for the evaporation. Such a "never-come-back" approach is used by some 
other exciton models, for instance, by the Multistage Preequilibrium Model (MPM) used in 
LAHET pjl] and by FLUKA [I32]. This approach was used in the CEM2k [17] version of the 
CEM and it allowed us to describe much better the p-|-A reactions measured at GSI in inverse 
kinematics at energies around 1 GeV/nucleon. Nevertheless, the "never-come-back" approach 
seems unphysical, therefore we no longer use it. We now address the problem of emitting fewer 
preequilibrium particles in the CEM by following Veselsky jl33] . We assume that the ratio of 
the number of quasi-particles (excitons) n at each preequilibrium reaction stage to the number 
of excitons in the equilibrium configuration neq, corresponding to the same excitation energy, 
to be a crucial parameter for determining the probability of preequilibrium emission 
probability for a given preequilibrium reaction stage is evaluated using the formula 



for n < Ueq and equal to zero for n > rieq. The basic assumption leading to Eq. (38) is that 
Ppre depends exclusively on the ratio n/rieq as can be deduced from the results of Bohning |134] 



f{E) = 1 - exp{-jE) . 



(37) 



a = 0.1463, (3 = -0.0716, and 7 = 0.0542 . 




pre 



(38) 
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where the density of particle-hole states is approximately described using a Gaussian centered 
at n^q. The parameter free parameter and we assume no dependence on excitation 

energy [133] . Our calculations of several reactions using different values of apre show that an 
overall reasonable agreement with available data can be obtained using apre = 0.4-0.5 (see Fig. 
11 in Ref. [IB])- In CEM03.0x, we choose the fixed value apre = 0.4 and use Eqs. (34,38) as 
criteria for the transition from the preequilibrium stage of reactions to evaporation, instead of 
using the "never-come-back" approach along with Eq. (34), as was done in CEM2k. 

Algorithms of many preequilibrium routines are changed and almost all these routines are 
rewritten, which has speeded up the code significantly. Finally, some bugs were fixed as previ- 
ously mentioned. 

Several examples with results from the preequilibrium model used in our current event 
generators are shown in Figs. 19-24. 

The energy spectra of secondary particles emitted from 96 MeV n + ^^Fe interactions 
presented in Fig. 17 show that the main contribution to the total spectra of complex particles 
as calculated by CEM03.01 comes from preequilibrium emission: The coalescence mechanism 
contributes only a few percent to the total spectra from these reactions, while evaporation is 
important only in the production of low energy particles, below ~ 30 MeV. Fig. 19 shows several 
more detailed spectra for the same reactions, namely, CEM03.01 results for double-differential 
spectra of d and t at eight different angles compared with the measured data [112] . Let us recall 
that CEM03.0x produces almost all these deuterons and tritons via preequilibrium emission. 
We see that CEM03.0x describes quite well these double-differential spectra, except the very 
high-energy tails of spectra at the most forward angles, where we should expect a contribution 
from direct processes like pick-up and knock-out, not considered so far in our models. 

Fig. 20 shows examples of angle-integrated energy spectra, energy-integrated angular dis- 
tributions, and double-differential spectra of nucleons and complex particles from a reaction 
induced by intermediate-energy protons, namely from 62.9 MeV p + Pb, as calculated by 
CEM03.01 and compared with the recent measurements by Guertin et al. |135j . CEM03.01 
produces complex particles from this reaction also mainly via preequilibrium emission, and we 
see that it describes these experimental spectra quite well too. 

Fig. 21 shows another example of double-differential spectra of complex particles from 
proton-induced reactions, at higher energies, namely CEM03.01 calculated "^He spectra at 20, 
40, 60, 100, 120, and 140 degrees from 160 MeV p + ^^Al, ^^Co, and ^^"^Au compared with the 
Cowley et al. data |136j . We see again a good agreement between our calculations and the 
measurements; the main contribution in the production of these "^He by CEM03.01 is again 
from preequilibrium emission. 

One more example of double-differential spectra of complex particles from proton-induced 
reactions, at higher energies, namely CEM03.01 calculated ^He and "^He spectra at 20, 90, and 
160 degrees from 210, 300, and 480 MeV proton-silver interactions compared with the Green 
and Korteling data [137J is presented in Fig. 22. Again a fairly good agreement between the 
calculations and the data can be seen, and again the main contribution in the production of 
these particles by CEM03.01 is from preequilibrium emission. 

Finally, Fig. 23 show an example of double-differential spectra of complex particles from 
neutron-induced reactions, at even higher energies, namely CEM03.01 calculated double-differential 
spectra of p, d, and t at 54, 68, 90, 121, and 164 degrees from interactions of 542 MeV neutrons 
with copper and bismuth compared with the measurements by Franz et al. |138] . We see again 
a good agreement between the calculations and the data and again the main contribution in the 
production of d and t from these reactions by CEM03.01 is from preequilibrium emission. Let 
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us note that we analyzed these data in 1992 [26], with the version of CEM we had at that time, 
CEM92m |78] . If we compare the agreement with the experimental data [138j of the results we 
got in 1992 [26] with our current CEM03.01 results, there is a tremendous improvement in the 
description of these data. 

To conclude this Section, Fig. 24 shows several examples of gas production cross-section 
calculations by CEM03.01 from n+^^Fe, p+^^^Pb, and n+^^^U reactions compared with experi- 
mental data [i39j-[llT] and results by TALYS from [HO], McGNASH [US], and IPPE-99 [113] . 
The overall agreement of the results by CEM03.01 with these data is similar to that achieved 
by TALYS, McGNASH, and IPPE-99. 

For the production of ^He from n+'^^^U (bottom plot in Fig. 24), we show not only the 
total gas production cross section, but also contributions to the total yield from preequilibrium 
emission, from events with evaporation that are not followed by fission, from evaporation before 
fission and from evaporation of fission fragments in events with fission, the total evaporation 
component (both with and without fission), and from coalescence, as calculated by CEM03.01. 
We see that the preequilibrium contribution to the total yield of ^He from this reaction is about 
one order of magnitude higher than contributions from other reaction mechanisms considered 
by our model. 




20 40 60 80 100 20 40 60 80 100 

T (MeV) T (MeV) 



Figure 19: Double-differential spectra of d and t at 20, 40, 60, 80, 100, 120, 140, and 160 
degrees from 96 MeV n + ^^Fe calculated by CEM03.01 compared with experimental data by 
Blideanu et al. [112] . 
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Figure 20: Experimental angle-integrated energy spectra (upper left plot), energy- integrated 
angular distributions (upper right plot), and double-differential spectra of nucleons and complex 
particles from 62.9 MeV p + Pb [135J compared with CEM03.01 results. 
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21: CEM03.01 calculated double-differential spectra of ^He at 20, 40, 60, 100, 120, and 
grees from 160 MeV p + ^''Al, ^^Co, and ^^''Au compared with the Cowley et al. data 
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Figure 22: CEM03.01 calculated double-differential spectra of ^He and "^He at 20, 90, and 160 
degrees from 210, 300, and 480 MeV proton-silver interactions compared with the Green and 
Korteling data [137j. 
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Figure 23: CEM03.01 calculated double-differential spectra of p, ci, and t at 54, 68, 90, 121, 
and 164 degrees from interactions of 542 MeV neutrons with copper and bismuth compared 
with the measurements by Franz et al. [138j . 
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Figure 24: Examples of H (p + and ^He total production cross sections calculated by 
CEM03.01 compared with experimental data |139j - p3T] and results by TALYS from [140j . 
McGNASH and IPPE-99 [US]. For the production of ^He from n+^^^U (bottom plot), 

we show the total cross section as predicted by CEM03.01 with a thick red solid line, as well 
as contributions to it from preequilibrium emission, from events with evaporation that are not 
followed by fission, from evaporation before fission, evaporation from fission fragments, total 
evaporation (both with and without fission), and from coalescence, shown with different color 
lines as indicated in the legend. 
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6. Evaporation 



CEM03.01 and LAQGSM03.01 and their later versions use an extension of the Generahzed 
Evaporation Model (GEM) code GEM2 by Furihata |144j - [T^B] after the preequilibrium stage 
of reactions to describe evaporation of nucleons, complex particles, and light fragments heavier 
than ^He (up to ^^Mg) from excited compound nuclei and to describe their fission, if the 
compound nuclei are heavy enough to fission {Z > 65). The GEM is an extension by Furihata 
of the Dostrovsky evaporation model jl47j as implemented in LAHET [131] to include up to 
66 types of particles and fragments that can be evaporated from an excited compound nucleus 
plus a modification of the version of Atchison's fission model |148j - [TSU] used in LAHET. Many 
of the parameters were adjusted by Furihata for a better description of fission reactions when 
using it in conjunction with the extended evaporation model. 

A very detailed description of the GEM, together with a large amount of results obtained 
for many reactions using the GEM coupled either with the Bertini or ISABEL INC models 
in LAHET may be found in [144^ I145j . Therefore, we present here only the main features 
of the GEM, following mainly [145] and using as well useful information obtained in private 
communications with Dr. Furihata. 

Furihata did not change in the GEM the general algorithms used in LAHET to simulate 
evaporation and fission. The decay widths of evaporated particles and fragments are estimated 
using the classical Weisskopf-Ewing statistical model [125] . In this approach, the decay prob- 
ability Pj for the emission of a particle j from a parent compound nucleus i with the total 
kinetic energy in the center-of-mass system between e and e + de is 

Pj{e)de = gj(Jinv{e) ^'^^^ — —ede, (39) 

where E [MeV] is the excitation energy of the parent nucleus i with mass Ai and charge Zi, 
and d denotes a daughter nucleus with mass and charge Z^ produced after the emission 
of ejectile j with mass Aj and charge Zj in its ground state, a^v is the cross section for 
the inverse reaction, pi and are the level densities [MeV]~^ of the parent and the daughter 
nucleus, respectively, gj = (2Sj + l)mj/TT^h'^, where 5*^ is the spin and rrij is the reduced 
mass of the emitted particle j. The Q-value is calculated using the excess mass M{A,Z) as 
Q = M{Aj, Zj) + M{Aii, Zd) — M[Ai, Zj). In GEM2, four mass tables are used to calculate Q 
values, according to the following priorities, where a lower priority table is only used outside 
the range of validity of the higher priority one: (1) the Audi-Wapstra mass table [151] . (2) 
theoretical masses calculated by Moller et al. [129J, (3) theoretical masses calculated by Comay 
et al. |152j . (4) the mass excess calculated using the old Cameron formula |153 ]. As does 
LAHET, GEM2 uses Dostrovsky's formula |147j to calculate the inverse cross section crj„^ for 
all emitted particles and fragments 

cyinv{e) = ffgO ( 1 + - ) , (40) 



e 



which is often written as 



<ygCn{l + h/e) for neutrons 

cTgCj^l — V/e) for charged particles 



where ag = nRl [fm^] is the geometrical cross section, and 

V = k.Z^Z^e^/R, (41) 
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is the Coulomb barrier in MeV. 

One important new ingredient in GEM2 in comparison with LAHET, which considers evap- 
oration of only 6 particles (n, p, d, t, ^He, and '^He), is that Furihata includes the possibility of 
evaporation of up to 66 types of particles and fragments and incorporates into GEM2 several 
alternative sets of parameters b, Cj, kj, R^, and Rc for each particle type. 

The 66 ejectiles considered by GEM2 for evaporation are selected to satisfy the following 
criteria: (1) isotopes with Zj < 12; (2) naturally existing isotopes or isotopes near the stability 
line; (3) isotopes with half-lives longer than 1 ms. All the 66 ejectiles considered by GEM2 are 
shown in Table 1. 



Table 1. The evaporated particles considered by GEM2 



z, 


Ejectiles 















n 














1 


P 


d 


t 










2 


3He 


^He 


^'He 


«He 








3 


6Li 


^Li 


8Li 


9Li 








4 


^Be 


^Be 


lOBe 


i^Be 


i^Be 






5 


8B 


lOg 


"B 


12b 


13b 






6 


IOq 




12q 


ISq 


14o 


15o 


160 


7 


12N 


13N 


14N 




16N 


17N 




8 


14Q 


150 


160 


1^0 


180 


190 


20O 


9 




18p 


19p 


20p 


21p 






10 


i^Ne 


i^Ne 


20Ne 


2iNe 


22Ne 


23Ne 


24Ne 


11 


2iNa 


22Na 


23Na 


24Na 


25Na 






12 


22Mg 


23Mg 


24Mg 


25Mg 


26Mg 


27Mg 


28Mg 



GEM2 includes several options for the parameter set in expressions (40,41): 

1) The "simple" parameter set is given as c„ = Cj = kj = 1, b = 0, and R^ = Rc = 
ro{AY^ + A^J^) [fm]; users need to input tq. 

2) The "precise" parameter set is used in GEM2 as the default, and we use this set in our 
present work. 

A) For all light ejectiles up to a [Aj < 4), the parameters determined by Dostrovsky et al. 
[W] are used in GEM2, namely: c„ = 0.76 + CaA"^'^ , b = (baA^^^^ - 0.050)/(0.76 + CaA^^^^) 
(and b = for Ad > 192), where Ca = 1.93 and ba = 1.66, Cp = 1 + c, Cd = I + c/2, q = 1 + c/3, 
C3jjc = Ca = 0, kp = k, kd = k + 0.06, kt = k + 0.12, kajj^ = ka — 0.06, where c, k, and ka 
are listed in Table 2 for a set of Zd- Between the Zd values listed in Table 2, c, k, and k^ are 
interpolated linearly. The nuclear distances are given by Rb = 1.5A^^^ for neutrons and protons, 
and 1.5(Ay^ + A]^^) for d, t, ^He, and a. 



Table 2. k, ka, and c parameters used in GEM2 



Zd 


k 


ka 


c 


< 20 


0.51 


0.81 


0.0 


30 


0.60 


0.85 


-0.06 


40 


0.66 


0.89 


-0.10 


> 50 


0.68 


0.93 


-0.10 
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The nuclear distance for the Coulomb barrier is expressed as Rc = Rd + Rj, where Rd = 
TqA^/^, Tq = 1.7, and Rj = for neutrons and protons, and Rj = 1.2 for d, t, ^He, and ^He. We 
note that several of these parameters are similar to the original values published by Dostrovsky 
et al. [147] but not exactly the same. Dostrovsky et al. |147j had = 2.2, ha = 2.12, and 
Tg = 1.5. Also, for the k, ka, and c parameters shown in Table 2, they had slightly different 
values, shown in Table 3. 



Table 3. kp, Cp, ka, and Ca parameters from Ref. |147] 



Zd 


kp 


Cp 


ka 


Ca 


10 


0.42 


0.50 


0.68 


0.10 


20 


0.58 


0.28 


0.82 


0.10 


30 


0.68 


0.20 


0.91 


0.10 


50 


0.77 


0.15 


0.97 


0.08 


> 70 


0.80 


0.10 


0.98 


0.06 



B) For fragments heavier than a {Aj > 4), the "precise" parameters of GEM2 use values 
by Matsuse et al. [154J, namely: cj = k = 1, Rb = Ro{Aj) + Ro{Ad) + 2.85 [fm], R^ = 
Ro{A^) + RoiAd) + 3.75 [fm], where Ro{A) = 1.12A^/^ - OMA-^/\ 

3) The code GEM2 contains two other options for the parameters of the inverse cross 
sections. 

A) A set of parameters due to Furihata for light ejectiles in combination with Matsuse's 
parameters for fragments heavier than a. Furihata and Nakamura determined kj for p, d, t, 
■^He, and a as follows [146] : 

kj = ci log(Zd) + C2 log{Ad) + C3. 
The coefficients ci, C2, and C3 for each ejectile are shown in Table 4. 

Table 4. ci, C2, and C3 for p, d, t, ^He, and a from |146j 



Ejectile 


Cl 


C2 


C3 


P 


0.0615 


0.0167 


0.3227 


d 


0.0556 


0.0135 


0.4067 


t 


0.0530 


0.0134 


0.4374 


3He 


0.0484 


0.0122 


0.4938 


a 


0.0468 


0.0122 


0.5120 



When these parameters are chosen in GEM2, the following nuclear radius R is used in the 
calculation of V and ag-. 

for A = 1 , 
for 2 < A < 4 , 
for 5 < A < 6 , 
for A = 7 , 
for A = 8 , 
for A = 9 , 
for A> 10 . 



R 





1.2 

2.02 

2.42 

2.83 

3.25 



iaua]/^ + 1 
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B) The second new option in GEM2 is to use Furihata's parameters for light ejectiles up to 
a and the Botvina et al. |155j parameterization for inverse cross sections for heavier ejectiles. 
Botvina et al. |155j found that ainv can be expressed as 

, _ , / (1 - ^/e) for e > + 1 [MeV], 

^in. ^9|exp[a(e- V- l)]/(\/ + l) for e < V + 1 [MeV], ^ ' 



where 



a = 0.869 + 9.91/Zj, 



V 



.^^2.173i±^i^^iii^[fm]. 

The expression of a^nv for e < \/ + 1 shows the fusion reaction in the sub-barrier region. 
When using Eq. (42) instead of Eq. (40), the total decay width for a fragment emission can 
not be calculated analytically. Therefore, the total decay width must be calculated numerically 
and takes much CPU time. 

The total decay width Tj is calculated by integrating Eq. (39) with respect to the total 
kinetic energy e from the Coulomb barrier V up to the maximum possible value, {E — Q). The 
good feature of Dostrovsky's approximation for the inverse cross sections, Eq. (40), is its simple 
energy dependence that allows the analytic integration of Eq. (39). By using Eq. (40) for (Tinv, 
the total decay width for the particle emission is 

r, = ^ r^^fl + -)p.iE - g - e)de. (43) 
Pi{E) Jv ^ 

The level density p{E) is calculated in GEM2 according to the Fermi-gas model using the 
expression [156] 

where a is the level-density parameter and 6 is the pairing energy in MeV. As does LAHET, 
GEM2 uses the 6 values evaluated by Cook et al. [157] . For those values not evaluated by 
Cook et al, 6's from Gilbert and Cameron |156] are used instead. The simplest option for the 
level-density parameter in GEM2 is a = Ad/8 [MeV~^], but the default is the Gilbert-Cameron- 
Cook-Ignatyuk (GCCI) parameterization from LAHET [131] : 

^ 1 - e"" / 1 - e""\ 

a = a ha/ 1 , (45) 

u \ ^ J 

where u = 0.05{E — 6), and 

aj = (0.1375 - 8.36 x 10"^^^) x A^, 

Ad/8 for Zd < 9 or Nd< 9, 

Ad{a' + 0.00917^) for others. 

For deformed nuclei with 54 < < 78, 86 < < 98, 86 < iV^ < 122, or 130 < Nd < 150, 
a' = 0.12 while a' = 0.142 for other nuclei. The shell corrections S is expressed as a sum of 
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separate contributions from neutrons and protons, i.e. S = S(Zd) + S{Nd) from [ISGj 1157] and 
are tabulated in |144j . 

The level density is calculated using Eq. (44) only for high excitation energies, E > E^, 
where E^ = + S and = 2.5 + 150/^^ (all energies are in MeV). At lower excitation 
energies, the following [156] is used for the level density: 

piE) = ^^eMiE-E,)/T), (46) 



where T is the nuclear temperature defined as 1/T = \/ ajV^ — l.'^/Ux- To provide a smooth 
connection of Eqs. (44) and (46) aX E = E^, Eq is defined as Eq = E^ — T(logT — 0.25 log a — 
1.251ogf/^ + 2v4f4). 

For E — Q — V < E^, substituting Eq. (46) into Eq. (44) we can calculate the integral 
analytically, if we neglect the dependence of the level-density parameter a on E: 

r, = t) + (f^ + ^)^o(t)}, (47) 

where /o('^) ^^nd /i(t, t^:) are expressed as 

m = e-^«/^(e*-l), 
hit,t,) = e-^«/^T{(t-t, + l)e*^-t-l}, 

where t= {E-Q-V)/T and = E^/T. For E-Q-V > E^, the integral of Eq. (43) cannot 
be solved analytically because of the denominator in Eq. (44). However, it is approximated as 

r, = ^r7^[A(t, Q + his, s,)e^ + (/3 + V){Io{Q - his, s.)e^}], (48) 

where h{s,Sx) and 13(5,53;) are given by 

his, s,) = 2V2{s~'/' + 1.5s-'/' + 3.75s-'/' - is''/' + 1.5s-'/' + 3.75s-' /')e''^-'}, 

his, s,) = iV2a)-'[2s-'/' + 4S-3/2 + 13.5s-'/' + GO.Os-'/' + 325.1255-^/2 

- iis' - sl)s-^/' + il.5s' + 0.5sl)s-'/' + (3.755^ + 0.255^)5^^/^ + (12.8755^ 

+ 0.625s^)s-^/2 ^ (59.06255^ + 0.93754)s;"/^ + (324.83^ + 3.28sl)s-^^/'}e'^-'], 

with s = 2^aiE - Q - V - 6) and = 2^a(E^ - 5). 

The particle type j to be evaporated is selected in GEM2 by the Monte-Carlo method 
according to the probability distribution calculated as Pj = Tj / Tj, where Tj is given by 
Eqs. (47) or (48). The total kinetic energy e of the emitted particle j and the recoil energy of 
the daughter nucleus is chosen according to the probability distribution given by Eq. (39). The 
angular distribution of ejectiles is simulated to be isotropic in the center-of-mass system. 

According to Friedman and Lynch |158j . it is important to include excited states in the 
particle emitted via the evaporation process along with evaporation of particles in their ground 
states, because it greatly enhances the yield of heavy particles. Taking this into consideration, 
GEM2 includes evaporation of complex particles and light fragments both in the ground states 
and excited states. An excited state of a fragment is included in calculations if its half-life 
Ti/2is) satisfies the following condition: 

ln2 r*' ^ ' 
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where T* is the decay width of the excited particle (resonance). GEM2 calculates T* in the 
same manner as for a ground-state particle emission. The Q-value for the resonance emission is 
expressed as Q* = Q + E*, where Ej is the excitation energy of the resonance. The spin state 
of the resonance S* is used in the calculation of gj, instead of the spin of the ground state Sj. 
GEM2 uses the ground state masses mj for excited states because the difference between the 
masses is negligible. 

Instead of treating a resonance as an independent particle, GEM2 simply enhances the 
decay width Tj of the ground state particle emission as follows: 

r. = r? + Er"' (so) 

n 

where r° is the decay width of the ground state particle emission, and Tj is that of the nth 
excited state of the particle j emission which satisfies Eq. (49). 

The total-kinetic-energy distribution of the excited particles is assumed to be the same 
as that of the ground-state particle. S*, E*, and T1/2 used in GEM2 are extracted from the 
Evaluated Nuclear Structure Data File (ENSDF) database maintained by the National Nuclear 
Data Center at Brookhaven National Laboratory [159] . 

Note that when including evaporation of up to 66 particles in GEM2, its running time 
increases significantly compared to the case when evaporating only 6 particles, up to ^He. The 
major particles emitted from an excited nucleus are n, p, d, t, ^He, and ^He. For most cases, 
the total emission probability of particles heavier than a is negligible compared to those for the 
emission of light ejectiles. Our detailed study of different reactions (see, e.g., [160] and references 
therein) shows that if we study only nucleon and complex-particle spectra or only spallation 
and fission products and are not interested in light fragments, we can consider evaporation 
of only 6 types of particles in GEM2 and save much time, getting results very close to the 
ones calculated with the more time consuming "66" option. In CEM03.01 and LAQGSM03.01, 
we have introduced an input parameter called nevtype that defines the number of types of 
particles to be considered at the evaporation stage. The index of each type of particle that can 
be evaporated corresponds to the particle arangement in Table 1, with values, e.g., of 1, 2, 3, 4, 
5, and 6 for n, p, d, t, '^He, and ^He, with succeeding values up to 66 for ^^Mg. All 66 particles 
that can possibly evaporate are listed in CEM03.01 and LAQGSM03.01 together with their 
mass number, charge, and spin values in the block data bdejc. For all ten examples of inputs 
and outputs of CEM03.01 included in Appendices 1 and 2 of the CEM03.01 User Manual [23] , 
whose results are plotted in the figures in Appendix 3 of [23], we have performed calculations 
taking into account only 6 types of evaporated particles (nevtype = 6) as well as with the "66" 
option (nevtype = 66) and we provide the corresponding computing time for these examples 
in the captions to the appropriate figures shown in Appendix 3 of Ref. [23]. The "6" option 
can be up to several times faster than the "66" option, providing meanwhile almost the same 
results. Therefore we recommend that users of CEM03.01 and LAQGSM03.01 use 66 for the 
value of the input parameter nevtype only when they are interested in all fragments heavier 
than ^He; otherwise, we recommend the default value of 6 for nevtype, saving computing time. 
Alternatively, users may choose intermediate values of nevtype, for example 9 if one wants to 
calculate the production of ^Li, or 14 for modeling the production of ^Be and lighter fragments 
and nucleons only, while still saving computing time compared to running the code with the 
maximum value of 66. 

Examples of calculation by CEM03.01 the reactions 800 MeV p + ^^"^Au and 1 GeV p + 
^^Fe using the "6" and "66" options are shown in Figs. 25 and 26. 
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Figure 25: The measured |161] mass and charge distributions of the product yields from the 
reaction 800 MeV/A ^^^Au+p and of the mean kinetic energy of these products, and the mass 
distributions of the cross sections for the production of thirteen elements with the charge Z 
from 20 to 80 (open symbols) compared with CEM03.01 results. The results shown in this 
figure are for ten million simulated inelastic events. The nevtype=66 option requires 53 hr 34 
min 31 sec of computing time on a SunBlade 100, 500 MHz computer, while the nevtype=6 
option requires only 12 hr 48 min 3 sec, providing almost the same results for the spallation and 
fission products. The fragment {2 < Z < 13, 6 < A < 29) results are very different, therefore 
we need to use the option nevtype=66 when we are interested in fragment production. 
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Figure 26: Experimental mass distributions of the yields of eight isotopes from Na to Mn 
[162] and of all light fragments from Li to O [163j from the reaction 1 GeV/A ^^Fe+p and the 
charge distribution of the product yield compared with CEM03.01 results. The results shown 
in this figure are for ten million simulated inelastic events. The nevtype=66 option requires 
27 hr 38 min 37 sec of computing time on a SunBlade 100, 500 MHz computer, while the 
nevtype=6 option requires only 8 hr 27 min 48 sec, providing almost the same results for the 
spallation products. The yields of light fragments, especially of Li and Be, differ by several 
orders of magnitude, therefore we need to use the option nevtype > 6 when we are interested 
in light-fragment production. 
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Let us note that increasing the number of types of particles that can be evaporated from 6, 
employed by all old versions of CEM and LAQGSM which use our own evaporation model based 
on the Weisskopf-Ewing statistical theory of particle evaporation [125j and consider some newer 
features according to Ref. |128] (see details, e.g., in [H]), to 66 in the 03.01 and later versions, 



where we replaced our evaporation model with a modification of GEM2 by Furihata |144] - p^5] 
as described above, together with considering the Fermi breakup model to disintegrate exited 
nuclei with A < 13 in the new versions of our codes, allows us to improve considerably the 
description of fragment yields from various nuclear reactions. Two examples with such results 
are shown in Figs. 27 and 28. 

Fig. 27 shows an example of an excitation function of interest to space applications used 
to estimate the risk assessment for electronic devices (SEUs and MBUs), namely ^^Si(p,x)^Be, 
calculated with an already quite old version of CEM, CEM2k [17], and with our recent versions 
CEM03.01 [23] and CEM03.02 |63] compared with experimental data from our T-16 compilation 
[164] . CEM2k, which does not consider evaporation of up to 66 types of particles and the Fermi 
breakup of light excited nuclei, underestimates these data by two orders of magnitude, providing 
only a tiny production of ^Be, only at high incident energies as final residual nuclei produced 
at the end of the reaction after INC, preequilibrium, and evaporation. The new versions 
of our codes describe this reaction much better than CEM2k, though there is still room for 
improvement and we plan a further improvement of the description of fragment production by 
our codes. 
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Figure 27: 
CEM03.01 



Comparison of the ^^Si(p,X)^Be excitation function calculated by CEM2k |17j . 
[23j . and CEM03.02 [63] with experimental data from our T-16 compilation [164] . 



Fig. 28 compares results from CEM03.01 [23j and from its "S" and "G" modifications 
described briefly in Section 10, CEM03.S1 and CEM03.G1 [97], for the total production cross 
sections of Li and Be isotopes produced in interactions of 1.2 GeV protons with thirteen target 
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nuclei from Al to Th, measured recently at the Cooler Synchrotron Facility COSY of the 
Forschungszentrum Jiilich [165] . 
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Figure 28: Comparison of measured |165 ] (symbols) production cross sections of Li and Be 
isotopes with kinetic energies below 100 MeV for 1.2 GeV proton-induced reactions on tar- 
gets between Al and Th with results from CEM03.01, CEM03.S1, and CEM03.G1 (lines), as 
indicated. 



We see that on the whole, with only a few exceptions, all versions of our current codes 
describe reasonably the shape and the absolute values of the measured total production cross 
sections for all these fragments, while the old versions would produce such fragments only as 
residual nuclei and would underestimate the data by orders of magnitude, as discussed above. 
More details on comparisons of our results with the NESSI data [165] may be found in Ref. 

m- 

However, a reasonable description by our current codes of the experimental integrated frag- 
ment yields demonstrated above does not yet assure a similar agreement with the fragment 
spectra, where we still have some serious problems to solve. Fig. 29 shows examples of some 
problems we still have in CEM03.03, CEM03.01 [23], as well as in the "S" and "G" versions 
CEM03.S1 and CEM03.G1 ^ in the description of ^Li spectra at 20, 45, 60, 90, and 110 
degrees from proton-aluminum interactions at 200 MeV measured recently in Ref. |166] (we 
got very similar results also for the ''Li, ''Be, ^Be, ^°Be, and ^°B spectra). We see that on the 
whole, with only a few exceptions, all these recent versions of our codes describe reasonably the 
low-energy part of the measured spectra, where fragments are produced mostly via evaporation. 
But all codes fail to reproduce correctly the high-energy tails of the spectra of all the fragments 
heavier than ^He. Note that the authors of this experiment jT66] have measured similar frag- 
ment spectra also from ^^Co and ^^^Au, and the problems with their correct description by 
our codes for these systems are only graver. We believe that we have these problems (just like 
other models have; as of today, we do not know any models or codes able to reproduce well all 
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these data) because the high-energy fragments from these reactions are produced via preequi- 
hbrium mechanisms of nuclear reactions, not considered by our current models for fragments 
heavier than ^He. Our models consider preequilibrium emission of particles only up to ^He, 
this is probably the reason why they describe quite well the ^He spectra, and not so well but 
still reasonably the ^He spectra, but fail completely to reproduce the high energy tails of all 
heavier fragments. The incident energy of bombarding protons in these reactions is too low 
to produce fragments via multigragmentation; this is the reason why the "S" version does not 
work well here. The "G" version describes production of fragments from these reactions only 
via fission-like binary decays; this mechanism contributes only to the evaporation part of the 
spectra, therefore it also fails to reproduce the high-energy results (see more details on the "G" 
and "S" versions of our codes in Section 10). We plan to address preequilibrium emission of 
fragments heavier than ^He in the future, with a hope to solve these problems. 

7. Fission 

The fission model used in GEM2 is based on Atchison's model |148i 1149] as implemented in 
LAHET |131j . often referred in the literature as the Rutherford Appleton Laboratory (RAL) 
fission model, which is where Atchison developed it. In GEM2 there are two choices of param- 
eters for the fission model: one of them is the original parameter set by Atchison |148[ 1149] as 
implemented in LAHET |131j . and the other is a parameter set developed by Furihata [144(1145] . 

7.1. Fission Probability 

The Atchison fission model is designed to describe only fission of nuclei with Z > 70. It 
assumes that fission competes only with neutron emission, i.e., from the widths Tj of n, p, 
d, t, ^He, and ^He, the RAL code calculates the probability of evaporation of any particle. 
When a charged particle is selected to be evaporated, no fission competition is taken into 
account. When a neutron is selected to be evaporated, the code does not actually simulate 
its evaporation, instead it considers that fission may compete, and chooses either fission or 
evaporation of a neutron according to the fission probability Pf. This quantity is treated by 
the RAL code differently for the elements above and below Z = 89. The reasons Atchison split 
the calculation of the fission probability Pf are: (1) there is very little experimental information 
on fission in the region Z = 85 to 88, (2) the marked rise in the fission barrier for nuclei with 
Z'^/A below about 34 (see Fig. 2 in [149] ) together with the disappearance of asymmetric mass 
splitting, indicates that a change in the character of the fission process occurs. If experimental 
information were available, a split between regions around Z'^/A ^ 34 would be more sensible 

1) 70 < Zj < 88. For fissioning nuclei with 70 < Zj < 88, GEM2 uses the original 
Atchison calculation of the neutron emission width and fission width Fj to estimate the 
fission probability as 

= iTTt; = i + F„,/F/ ^^^^ 

Atchison uses |148l 1149] the Weisskopf and Ewing statistical model |125] with an energy- 
independent pre-exponential factor for the level density (see Eq. (44)) and Dostrovsky's |147] 
inverse cross section for neutrons and estimates the neutron width F„ as 
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T (MeV) 

Figure 29: Experimental ^He and ^Li spectra at 20, 45, 60, 90, and 110 degrees from proton- 
aluminum interactions at 200 MeV [l66j compared with results from CEM03.03, CEM03.01, 
CEM03.S1, and CEM03.G1, as indicated. 



54 



r„ = 0.352(1.68Jo + l.93Ay^Ji + Af^^{0.76Ji - O.OSJq)), (52) 
where Jq and Ji are functions of the level-density parameter a„ and s„(= 2^ an{E — Qn — 6)), 

(sn - l)e^" + 1 



Jo 

Ji = - 



2a„ 

(24 - Qsn + 6)e^" + 4-6 



Note that the RAL model uses a fixed value for the level-density parameter a^, namely 

a„ = (A, - l)/8, (53) 

and this approximation is kept in GEM2 when calculating the fission probability according to 
Eq. (51), although it differs from the GCCI parameterization (45) used in GEM2 to calculate 
particle evaporation widths. The fission width for nuclei with 70 < Zj < 88 is calculated in the 
RAL model and in the GEM as 

r,.^Ii^}}£L±l, ,54) 



where sj = 2^/(2/ (-E — -B/ — 6) and the level-density parameter in the fission mode a/ is fitted 
by Atchison to describe the measured Tf/Tn to be [149] : 

aj = a„ (^1.08926 + 0.01098(x - 31.08551)^) , (55) 

and X = Z^/A. The fission barriers Bf [MeV] are approximated by 

Z^ /^2\2 

5/ = Q„ + 321.2-16.7-^ + 0.218f -^j . (56) 

Note that neither the angular momentum nor the excitation energy of the nucleus are taken 
into account in finding the fission barriers. 

2) Zj > 89. For heavy fissioning nuclei with Zj > 89 GEM2 follows the RAL model plSllTig] 
and does not calculate at all the fission width Fj and does not use Eq. (51) to estimate the 
fission probability Pf. Instead, the following semi-empirical expression obtained by Atchison 
|148t I149j by approximating the experimental values of F„/Fj published by Vandenbosch and 
Huizenga |167] is used to calculate the fission probability: 

log(F„/Fy) = C{Z,){A, - Ao(Z,)), (57) 

where C(Z) and Ao{Z) are constants depending on the nuclear charge Z only. The values of 
these constants are those used in the current version of LAHET |131j and are tabulated in 
Table 5 (note that some adjustments of these values have been done since Atchison's papers 
[UHlIIin] were pubfished). 

In this approach the fission probability Pf is independent of the excitation energy of the 
fissioning nucleus and its angular momentum. 
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1.2. Mass Distribution 



The selection of the mass of the fission fragments depends on whether the fission is symmetric 
or asymmetric. For a pre- fission nucleus with Zf /Ai < 35, only symmetric fission is allowed. For 
Zf/Ai > 35, both symmetric and asymmetric fission are allowed, depending on the excitation 
energy of the fissioning nucleus. No new parameters were determined for asymmetric fission in 
GEM2. 

For nuclei with Zf/Ai > 35, whether the fission is symmetric or not is determined by the 
asymmetric fission probability Pasy 

P-^y = i + 4870e-o-36iJ- 
Table 5. C{Z) and Ao{Z) values used in GEM2 
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0.17123 


275.30 


106 


0.17123 


278.80 



7. 2. a. Asymmetric fission. For asymmetric fission, the mass of one of the post-fission fragments 
Ai is selected from a Gaussian distribution of mean A^ = 140 and width Gm = 6.5. The mass 
of the second fragment is A2 = A^ — Ai. 

7.2.b. Symmetric fission. For symmetric fission, Ai is selected from a Gaussian distribution of 

mean Aj = Ai/2 and two options for the width ctm as described below. 
The first option for choosing ctm is the original Atchison approximation: 

_ (3.97 + A25{E-Bf)- 0m212{E-Bfy, 
^^"\ 25.27, 

for (E — Bf) below or above 100 MeV, respectively. In this expression all values are in MeV 
and the fission barriers Bf are calculated according to Eq. (56) for nuclei with Zi < 88. For 
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nuclei with Zi > 88, the expression by Neuzil and Fairhall |168] is used: 



Bf = C-OM{Zf/A,), (60) 

where C = 18.8,18.1,18.1, and 18.5 [MeV] for odd-odd, even-odd, odd-even, and even-even 
nuclei, respectively. 

The second option in GEM2 for (Tm (used here) was found by Furihata [1441 1145] as: 

aM = C3{Z^/A,y + C^{Zf/A,) + C^iE - Bf) + C,. (61) 

The constants Cg = 0.122, C4 = -7.77, C5 = 3.32 x lO^^, and Cq = 134.0 were obtained by 
fitting with GEM2 the recent Russian collection of experimental fission-fragment mass distri- 
butions [169j . In this expression, the fission barriers Bf hj Myers and Swiatecki |170] are used. 
More details may be found in Ref. [145] . 

7.3. Charge Distribution 



The charge distribution of fission fragments is assumed to be a Gaussian distribution of 
mean Zf and width az- Zf is expressed as 



Zi + Z[ — Z'2 



(62) 



where 



Z[ 



65.5A 
131 + 



/ = 1 or 2. 



(63) 



The original Atchison model uses az = 2.0. An investigation by Furihata [145] suggests that 
az = 0.75 provides a better agreement with data; therefore az = 0.75 is used in GEM2 and in 
our code. 



7.4- Kinetic Energy Distribution 

The kinetic energy of fission fragments [MeV] is determined by a Gaussian distribution with 
mean e/ and width a^^. 

The original parameters in the Atchison model are: 

€f = 0.133Z,VA^^^ - 11-4, 
a,^. = 0.084e/. 

Furihata's parameters in the GEM, which we also use, are: 

0.131Zf/Ay\ 
O.lOiZf/Ay^ + 24.3, 



for Zf/A^^ < 900 and 900 < Zf/Ay^ < 1800, respectively, according to Rusanov et al. [ISS]. 
By fitting the experimental data by Itkis et al. [171] . Furihata found the following expression 
for a^j 

,^^, = |c;i(m^/^- 1000)+C„ 
57 



for ZfjA^ above and below 1000, respectively, and the values of the fitted constants are 
C\ = 5.70 X 10~^ and C2 = 86.5. The experimental data used by Furihata for fitting are the 
values extrapolated to the nuclear temperature 1.5 MeV by Itkis et al. p7lj . More details may 
be found in [145j . 

We note that Atchison has also modified his original version using recent data and published 
1 172] improved (and more complicated) parameterizations for many quantities and distributions 
in his model, but these modifications [T72| have not been included either in LAHET or in GEM2. 

7.5. Modifications to GEM2 m CEM03.01 and LAQGSM03.01 

First, we fixed several observed uncertainties and small errors in the 2002 version of GEM2 
Dr. Furihata kindly sent us. Then, we extended GEM2 to describe fission of lighter nuclei, 
down to Z > 65, and modified it [20] so that it provides a good description of fission cross 
sections when it is used after our INC and preequilibrium models. 

If we had merged GEM2 with the INC and preequilibrium-decay modules of GEM or of 
LAQGSM without any modifications, the new code would not describe correctly fission cross 
sections (and the yields of fission fragments). This is because Atchison fitted the parameters of 
his RAL fission model when it followed the Bertini ING ^173] which differs from ours. In addi- 
tion, Atchison did not model preequilibrium emission. Therefore, the distributions of fissioning 
nuclei in A, Z, and excitation energy E* simulated by Atchison differ significantly from the 
distributions we get; as a consequence, all the fission characteristics are also different. Furihata 
used GEM2 coupled either with the Bertini ING [173] or with the ISABEL [Hi] ING code, 
which also differs from our ING, and did not include preequilibrium particle emission. Therefore 
the distributions of fissioning nuclei simulated by Furihata differ from those in our simulations, 
so the parameters adjusted by Furihata to work well with her ING are not appropriate for us. 
To get a good description of fission cross sections (and fission-fragment yields) we have modified 
at least two parameters in GEM2 as used in GEM03.01 and LAQGSM03.01 (see more details 

in onHisg). 

The main parameters that determine the fission cross sections calculated by GEM2 are the 
level-density parameter in the fission channel, aj (or more exactly, the ratio aj /an as calculated 
by Eq. (55)) for preactinides, and parameter C{Z) in Eq. (57) for actinides. The sensitivity 
of results to these parameters is much higher than to either the fission-barrier heights used 
in a calculation or other parameters of the model. Therefore we choose [20j to adjust only 
these two parameters in our merged codes. We do not change the form of systematics (55) 
and (57) derived by Atchison. We only introduce additional coefficients both to a/ and C(Z), 
replacing a/ — Ca x a/ in Eq. (55) and C{Z.i) ^ Cc y< C{Zi) in Eq. (57) and fit Ca and Cc to 
experimental proton-induced fission cross sections covered by Prokofiev's systematics [175] for 
both GEM03.01 and LAQGSM03.01. No other parameters in GEM2 have been changed. For 
preactinides, we fit only Ca- The values of Ca found in our fit to Prokofiev's systematics are 
close to one and vary smoothly with the proton energy and the charge or mass number of the 
target. This result gives us some confidence in our procedure, and allows us to interpolate the 
values of Ca for nuclei and incident proton energies not analyzed by Prokofiev. For actinides, 
as described in [13, US], we have to fit both Ca and Cc- The values of Ca we find are also very 
close to one, while the values of Cc are more varied, but both of them change smoothly with 
the proton energy and Z or A of the target, which again allows us to interpolate them for nuclei 
and energies outside Prokofiev's systematics. 

We fix the fitted values of Ca and Cc in data blocks in our codes and use the routines fitafpa 
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and fitafac to interpolate to nuclei not covered by Prokofiev's systematics. We believe that 
such a procedure provides a reasonably accurate fission cross-section calculation, at least for 
proton energies and target nuclei not too far from the ones covered by the systematics. 

Let us mention that the situation with the fitting procedure of parameters Ca and Cc is 
quite tricky, as it should be redone after all major improvements of other parts of our codes 
describing INC, preequilibrium, or evaporation. This is a major minus of such types of models 
like GEM2 that are based mostly on systematics of available experimental data rather than 
on fundamental physics. (This was the main reason we started to develop our own improved 
evaporation and fission models for our codes [12] that would describe experimental data not 
worse than GEM2 but would be based on physics rather than on systematics of available data; 
this work is not completed yet, so we have to use GEM2 in our codes so far.) Indeed, after 
making our major improvements to the INC and preequilibrium parts CEM and LAQGSM as 
described above, the mean values of the mass and charge numbers, A and Z of the excited 
compound nuclei produced after the preequilibrium stage of nuclear reactions and their mean 
excitation energy E* have changed slightly, which affects the probability of heavy compound 
nuclei (especially preactinides) to fission. This means that the procedure of fitting the Ca 
and Cc parameters which we performed in [20| to provide the best description by CEM2k and 
LAQGSM of fission cross sections was no longer correct. We had to redo this fitting for the 
latest versions of our CEM03.01 and LAQGSM03.01, ensuring that they describe as well as 
possible fission cross sections from various reactions. Fig. 30 shows examples of fission cross 
sections for proton-induced reactions on ^^^W, ^^^W, ^^^W and ^^^W. The improved CEM03.01 
reproduces the recent Uppsala measurements |176 j of proton-induced fission cross sections. 
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Figure 30: Experimental |176j proton-induced fission cross sections of ^^^W, ^^^W, ^'^^W, and 
^^^W compared with improved (red solid lines) and old (blue dashed lines, from |176j ) CEM03.01 
calculations. 

Another example of some problems of the fission model of GEM2 based on systematics 
of measured data is shown in Fig. 31 for fragment mass distributions from several reactions 
induced by bremsstrahlung with the maximum energy Eq of 15, 20, 30, and 70 MeV on ^ss-q ^^-^^ 
^^^U targets measued in |177] compared with results by CEM03.02 [63]. We see that the fission 
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model of GEM2 used in our CEM03.02 describes surprisingly well these mass distributions of 
fission fragments at Eq = 15, 20, and 30 MeV, where experimentasl data were available to fit the 
mass distribution of fission fragments in the GEM2 systematics, but fails to reproduce correctly 
such distribution at £"0 = 70 MeV, requiring an additional refitting of the GEM systematics at 
this energy to improve the agreement with the data. 




70 80 90 100 110 120 130 140 150 160 170 
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Figure 31: Experimental |177] fission fragment mass distributions (symbols) from reactions 
induced by bremsstrahlung with the maximum energy Eq of 15, 20, 30, and 70 MeV on ^ss-q 
and 238-[j targets compared with results by CEM03.02 [63] (histogams). 



Finally, to conclude the fission Section, Figs. 32 and 33 shows examples of spallation, fission, 
and fragmentation experimental product cross sections from 1 GeV p + ^38^ ji07[ HQgj and 2 
GeV d + ^°^Pb |178j measured recently at GSl in inverse kinematics compared with results by 
CEM03.01 {p + U, Fig. 32) and LAQGSM03.01 {d + Pb, Fig. 33), respectively. Our models 
describe quite well the fission- fragment yields, as well as the spallation and fragmentation 
product cross sections and agree with most of the GSI data with an acuracy of a factor of two 
or better. Similar agreements are obtained with other reactions measured at GSI in inverse 
kinematics (see, e.g., [311 H2l IT601 fT79] ) . 
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Figure 32: Comparison of measured [1071 I108j spallation, fission, and fragmentation cross 
sections for the reaction ^^^U(l GeV/A) + p (filled circles) with our CEM03.01 results (open 
circles). Experimental data for isotopes from B to Co, from Tb to Ta, and for Np and Pu are 
not available so we present only our predictions. 
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Figure 33: Comparison of measured [178] spallation and fission-product cross sections for the 
reaction 208pb(l GeV/A) + d (filled circles) with our LAQGSM03.01 results (open circles). 



8. The Fermi Breakup Model 

After calculating the coalescence stage of a reaction, CEM03.01 and LAQGSM03.01 move 
to the description of the last slow stages of the interaction, namely to preequilibrium decay and 
evaporation, with a possible competition of fission. But as mentioned above, if the residual 
nuclei have atomic numbers with A < 13, CEM03.01 and LAQGSM03.01 use the Fermi breakup 
model |18Uj to calculate their further disintegration instead of using the preequilibrium and 
evaporation models. The newer 03.02 versions of our codes use the Fermi breakup model also 
during the preequilibrium and/or evaporation stages of reactions, when the residual nucleus 
has an atomic number with A < 13. 

Finally, the latest 03.03 versions of our codes use the Fermi breakup model also to disinte- 
grate the unstable fission fragments with A < 13 that can be produced in very rare cases of 
very asymmetric fission. 

All formulas and details of the algorithms used in the version of the Fermi breakup model 
developed in the former group of Prof. Barashenkov at Joint Institute for Nuclear Research 
(JINR), Dubna, Russia and used by our codes may be found in [53]. All the information 
needed to calculate the breakup of an excited nucleus is its excitation energy U and the mass 
and charge numbers A and Z. The total energy of the nucleus in the rest frame will be 
E = U + M{A, Z), where M is the mass of the nucleus. The total probability per unit time 
for a nucleus to break up into n components in the final state {e.g., a possible residual nucleus, 
nucleons, deuterons, tritons, alphas, etc.) is given by 

W{E,n) = {V/nr-'pn{E), (66) 
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where p„ is the density of final states, V is the volume of the decaying system and Q = (27r^)^ 
is the normalization volume. The density Pn{E) can be defined as a product of three factors: 

Pn{E) = Mn{E)SnGn. (67) 

The first one is the phase space factor defined as 

/+00 p+oo / " \ / " / \ " 

■■■/ 6[J2pt]s[E-J2^P' + <]ll^'Pb, (68) 
-~ \b=l J \ 6=1 / 6=1 

where pb are fragment momenta. The second one is the spin factor 

n 

Sn = l[{2sb + l), (69) 
b=i 

which gives the number of states with different spin orientations. The last one is the permuta- 
tion factor 

k 

Gn = ll:^, (70) 
j=i r 

which takes into account identical particles in the final state {rij is the number of components 
of j'-type particles and k is defined by n = Yl^=i '^j)- -^^^ example, if we have in the final state 
six particles (n = 6) and two of them are alphas, three are nucleons, and one is a deuteron, 
then = l/(2!3!l!) = 1/12. For the non-relativistic case, the integration in Eq. (68) can be 
evaluated analytically (see, e.g., [53]) and the probability for a nucleus to disintegrate into n 
fragments with masses rrib, where 5 = 1, 2, 3, . . . , n is 



W{E, n) = SA (^) (^^^ n ' (71) 
V^y \Y.b=i^bl^^ ) r(3(n-l)/2) 



where r(x) is the gamma function. 

The angular distribution of n emitted fragments is assumed to be isotropic in the cm. 
system of the disintegrating nucleus and their kinetic energies are calculated from momentum- 
energy conservation. The Monte-Carlo method is used to randomly select the decay channel 
according to probabilities defined by Eq. (71). Then, for a given channel, the code calculates 
kinematic quantities for each fragment according to the n-body phase space distribution using 
Kopylov's method [181] . Generally, the Fermi breakup model considers formation of fragments 
only in their ground and those low-lying states which are stable for nucleon emission. However, 
as already mentioned in Section 2, several unstable fragments with large lifetimes: ^He, ^Li, 
®Be, ^B, etc. are considered as well by the initial version of the Fermi breakup model code 
as described in Ref. The randomly chosen channel will be allowed to decay only if the 

total kinetic energy £^fcj„ of all fragments at the moment of breakup is positive, otherwise a new 
simulation will be performed and a new channel will be selected. The total kinetic energy E^i^ 
can be calculated according to the equation: 

n 

Ekin = U + M{A, Z) - Ecoulomb - J^i^b + Cfo), (72) 

6=1 
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where rrib and eb are masses and excitation energies of the fragments, respectively, and Ecouiomb 
is the Coulomb barrier for the given channel. It is approximated by 



where Ai, and are the mass number and the charge of the 6-th particle of a given channel, 
respectively. Vq is the volume of the system corresponding to normal nuclear density and 
V = kVo is the decaying system volume (we assume = 1 in our codes). 

Thus, the Fermi breakup model we use has only one free parameter, or Vq, the volume 
of the decaying system, which is estimated as follows: 



where we use tq = 1.4 fm. 

There is no limitation on the number n of fragments a nucleus may break up into in our 
version of the breakup model, in contrast to implementations in other codes, such as n < 3 in 
MCNPX, orn < 7 in LAHET. 

In comparison with its initial version as described in [53], the Fermi breakup model used 
in CEM03.02 and LAQGSM03.02 has been modified [63] to decay the unstable light fragments 
that were produced by the original code. As mentioned above, the initial routines that de- 
scribe the Fermi breakup model were written more than twenty years ago in the group of Prof. 
Barashenkov at JINR, Dubna, and unfortunately had some problems. First, those routines 
allowed in rare cases production of some light unstable fragments like ^He, ^Li, ^Be, ^B, etc. 
as a result of a breakup of some light excited nuclei. Second, they allowed very rarely even 
production of "neutron stars" (or "proton stars"), i.e., residual "nuclei" produced via Fermi 
breakup that consist of only neutrons (or only protons). Lastly, those routines could even crash 
the code, due to cases of division by 0. All these problems of the Fermi breakup model routines 
were addressed and solved in CEM03.02 [63]; the changes were then put in LAQGSM03.02 
[63] . Several bugs are also fixed. However, even after solving these problems and after im- 
plementing the improved Fermi breakup model into CEM03.02 and LAQGSM03.02 [63j, our 
event generators still could produce some unstable products via very asymmetric fission, when 
the excitation energies of those fragments were below 3 MeV so they were not checked and 
disintegrated with the Fermi breakup model. Our analysis [12] had shown that such events 
could occur very rarely, in less than 0.0006% of all simulated events, so that production of such 
unstable nuclides affects by less than 0.0006% the other correct cross sections calculated by 
our codes. However, these unstable nuclides are not physical and should be eliminated. This 
was the reason why a universal checking of all unstable light products has been incorporated 
into CEM03.03 and LAQGSM03.03. Such unstable products are forced to disintegrate via 
Fermi breakup independently of their excitation energy. The latest versions of the CEM03 and 
LAQGSM03 event generators do not produce any such unstable products. 

Examples of several results from CEM03.02 compared with results from CEM03.01 and 
from the older versions CEM95 and CEM2k and with experimental data for a number of 
reactions relevant to the Fermi breakup model are presented Figs. 34 to 36. Fig. 34 presents 
the mass distribution of the product yields from the reaction 730 MeV p -|- ^^Al calculated 
with CEM03.01 without considering the Fermi breakup mode during the preequilibrium and 
evaporation stages of reactions and with CEM03.02 that does consider this mode during these 
stages and uses an extended version of the model without unstable products, as described 




(73) 



V = 47ri?V3 = 47rr^A/3, 



(74) 
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above. CEM03.01 provides a small yield of unphysical unstable ^He (just like the older versions 
CEM95 and CEM2k do), while the newer version does not produce such unstable nuclides. The 
results from CEM03.02 also agree a little better with available experimental data (red squares 
on Fig. 34) than those from CEM03.01. Similar results are obtained for several other reactions 
on other light targets. 
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Figure 34: Mass distribution of the product yields from the reaction 730 MeV p + ^^Al calcu- 
lated with CEM03.01 without considering the Fermi breakup mode during the preequilibrium 
and evaporation stages of reactions (solid circles connected with a solid line) and with the im- 
proved version of the code CEM03.02 that does consider the Fermi breakup mode during the 
preequilibrium and evaporation stages of reactions (open circles connected with a dashed line), 
compared with experimental data available at a nearby energy of 800 MeV from the T-16 Lib 
compilation |164j (open red squares). 

Our results show that CEM03.02 and LAQGSM03.02 do not predict unstable unphysical 
nuclides and describe the yields of most products a little better than the 03.01 or older versions 
of GEM and LAQGSM. The question about how the spectra of different particles from different 
nuclear reactions are described by the 03.02 version in comparison with older versions of our 
codes is not obvious. We studied this question on several reactions and Fig. 35 presents examples 
of particle spectra calculated with GEM03.01 and GEM03.02 for the reaction 730 MeV p + Al. 
We see that spectra of n, p, d, t, '^He, and "^He calculated by GEM03.02 that uses the Fermi 
breakup model to describe disintegration of exited nuclei with A < 13 are almost the same as 
such spectra calculated by GEM03.01 which uses the preequilibrium and evaporation models 
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to describe such reactions when A > 12 after the INC stage of reactions. This result is very 
interesting from a physical point of view; it is not trivial at all and could not be forecast easily 
in advance. We see that spectra of these particles predicted by different models are almost the 
same. In other words, the theoretical spectra of particles do not depend much on the models we 
use to calculate them, but depend mostly on the final phase space calculated by these models: 
If the phase spaces calculated by different models are correct and near to each other, than the 
spectra of secondary particles calculated by these models would be also near to each other and 
would be not very sensitive to the dynamics of reactions considered (or not) by the models we 
use. 
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Figure 35: Angle-integrated energy spectra of secondary n, p, d, t, ^He, and '^He from the 
reaction 730 MeV p + ^^Al calculated with CEM03.01 without considering the Fermi breakup 
mode during the preequilibrium and evaporation stages of reactions (solid histograms) and 
with CEM03.02 that does consider the Fermi breakup mode during the preequilibrium and 
evaporation stages of reactions (dashed histograms). 



Fig. 36 shows examples of double differential experimental neutron spectra at 15, 30, 60, 
90, 120, and 150 deg from 3.0, 1.5, 0.8, and 0.597 GeV p + Al (symbols) compared with 
results by CEM03.01 without considering the Fermi breakup mode during the preequilibrium 
and evaporation stages of reactions (dashed histograms), with CEM03.02 that does consider 
this mode during the preequilibrium and evaporation stages of reactions (solid histograms), 
and with the old version of the code, CEM95 [9J. We see that the results from CEM03.02 
and CEM03.01 agree much better with the data at neutron energies of 20-100 MeV (where 
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the contribution from the preequihbrium mode is the most important) than the results from 
CEM95 do. Results from CEM03.02 for these neutron spectra are nearly identical to those 
from CEM03.01, just as observed for the energy spectra of n, p, d, t, ^He, and ^He shown in 
Fig. 35. 
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Figure 36: Experimental |182] - p^^ double differential neutron spectra at 15, 30, 60, 90, 
120, and 150 deg from 3.0, 1.5, 0.8, and 0.597 GeV p + Al (symbols) compared with results 
from CEM03.01 without considering the Fermi breakup mode during the preequihbrium and 
evaporation stages of reactions (dashed histograms), with CEM03.02 which does consider this 
mode during the preequihbrium and evaporation stages of reactions (solid histograms), and 
with the old version of the code, CEM95 [H]. 



9. Reactions Involving Pions and Photons; Excitation Functions 



Protons from the beam and energetic secondary particles generate in ADS targets and 
surrounding shielding not only nucleons, complex particles, and fragments heavier than "^He 
(in addition to residual nuclei), but also pions of different energies. The transport codes used 
in ADS applications should be able to describe correctly their spectra and yields, to transport 
them, and to consider their further interactions with nuclei. We have shown that our event 
generators describe reasonably well spectra and yields of nucleon and complex particles from 
various reactions, and a little worse but still reasonably the yields and low energy part of spectra 



67 



of fragments heavier than ^He, with still retaining some unsolved problems for the high energy 
part of fragment spectra. Here, we demonstrate that our codes describe reasonably well also 
pion spectra from different reactions as well as various pion-induced reactions. 

Fig. 37 shows an example of vr"*" spectra from 562.5 MeV n + Cu calculated by CEM03.01 
compared with the LANL measurements by Brooks et. al. [185J. We see a reasonably good 
agreement of CEM03.01 calculations with these pion spectra. Similar results are obtained 
with our codes for other pion spectra (and for the integrated pion yield) measured in various 
reactions induced by different projectiles at different energies (see, e.g., [22l[2llll2] and Fig. 12 
above). 
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Experimental 7r+ spectra from 562.5 MeV n + Cu [185] compared with CEM03.01 



Fig. 38 shows an example of neutron spectra from 1.5 GeV vr"*" + Fe calculated by CEM03.01 
compared with experimental data by Nakamoto et. al. [186J. We see again a reasonably good 
agreement of our results with these measured data. We obtain similar results for other pion- 
induced reactions. 

The secondary neutral pions produced in ADS targets and shielding materials would decay 
later into two photons, as their mean life is of only (8.4 ± 0.6) x 10"^^ s. That is, the trans- 
port codes should transport the produced photons and our event generators should be able to 
describe properly their further interactions with various nuclei. 

If the photon energy is above several GeV, such photonuclear interactions are calculated 
in transport codes using our high-energy event generator LAQGSM03.01. Several examples of 
results by LAQGSM03.01 for photonuclear reactions are shown in Fig. 14 of Section 3.2; more 
results for other high-energy photonuclear reactions may be found in Refs. [2H [88] . 

Interaction of photons of lower energies with nuclei are calculated in transport codes using 
our low-energy event generator CEM03.01. Figs. 39 and 40 show examples of proton spectra 
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Figure 38: Experimental neutron spectra from 1.5 GeV 7r+ + Fe |186j compared with CEM03.01 
results [23]. The results shown in this figure are for one million simulated inelastic events. The 
nevtype=66 option (see Section 6) requires 3 hr 30 min 52 sec of computing time on a 
SunBlade 100, 500 MHz computer, while the nevtype=6 option requires only 1 hr 51 min 58 
sec, providing almost the same results. 

from 300 MeV 7+^^Cu and isotopic yields of products produced by bremsstrahlung reactions 
on ^^^Au and ^°^Bi at Eq = 1 GeV calculated by CEM03.01 compared with experimental data 
|187[ 1188] . Figs. 39 and 40 demonstrate that CEM03.01 allows us to describe reasonably well 
many photonuclear reactions needed for ADS and other applications, as well as to analyze 
mechanisms of photonuclear reactions for fundamental studies (see [22j for more details). 

It is well known that the most difficult characteristics of nuclear reactions to be described 
by any theoretical model are excitation functions, i.e., cross sections for the production of 
a given isotope as functions of the energy of the projectile. We have studied with different 
versions of our models several thousands of excitations functions for various reactions (see, e.g., 
[38| |39| [76l 11891 11901 1191 j and references therein); therefore, we limit ourself here to analysis of 
only a few recent measurements published only recently |192j and calculated with CEM03.03 
just the day before our present Advanced Workshop. 

Figs. 41 to 43 show excitation functions for the production of ^He, ^He, ^^Ne, and ^^Ar from 
interactions of protons with Iron and Nickel measured recently by Ammon et al. |192j compared 
with our CEM03.03 resuhs, with results by TALYS tl93j and by INCL/ABLA [72l [T95] from 
[192] . with systematics by Konobeev and Korovin |194] (for '^He and ^He from p+Fe), with 
previous measurements by other authors, and with our predictions for the p+Fe excitation 
functions we published in 1997 [76]. On the whole, CEM03.03 is able to reproduce quite well 
these new measurements. 

It is interesting to observe on the right panels of Figs. 41 and 42 that the CEM95 [9J version 
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Figure 39: Proton spectra at 45°, 90°, and 135° from the reaction 300 MeV 7 + Cu. Symbols 
are experimental data from [187] and histograms are CEM03.01 results [23]. 

of CEM predicted these p+Fe excitation functions reasonably well: the old CEM95 predictions 
[76] agree quite well with the recently measured data [192] . 

We observe some disagreements of our results with the new data, especially for the produc- 
tion of ^^Ne at low energies from both Fe and Ni, indicating us that there are still problems to 
be solved in CEM03.03 and we have to improve further our models to describe properly this 
type of reaction. We see that TALYS [193] and INCL/ABLA ^|195] also encounter difficulties 
in reproducing these data, indicating similar problems occur with other models. 

10. CEM03.S1, CEM03.G1, LAQGSM03.S1, and LAQGSM03.G1 

From the results presented already and in the cited references, we can conclude that 
CEM03.01 and LAQGSM03.01 and their later 03.02 and 03.03 versions are able to predict well 
a large variety of nuclear reactions of interest to ADS. Therefore, they can be employed with 
confidence as reliable event generators in transport codes used as "workhorses" for ADS and 
other applications. This is true only if we are not interested in the production of intermediate- 
mass fragments from not too heavy targets at not too high energies, targets that do not fission 
from an "orthodox" point of view, therefore not producing such fragments. To be able to pre- 
dict production of intermediate-mass fragments from not too heavy targets at such intermediate 
energies with our codes, we still need to solve some problems, just as do authors of other similar 
Monte-Carlo codes. 

Such a problem is presented in the low-energy parts of the Fe(p,x)^^Ne and Ni(p,x)^^Ne 
excitation functions shown in Figs. 42 and 43; two more examples, for other reactions, are 
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Figure 40: Comparison of CEM03.01 results (green symbols) for the isotopic yields of products 
produced by bremsstrahlung reactions on ^^''Au and ^°^Bi at Eq = 1 GeV with experimental 
data (red symbols) from the review |188j and calculations by PICA3/GEM (blue symbols); the 
PICA3/GEM results are from several publications and are presented in Fig. 18 of [188] with 
the corresponding citations. The mass yields for the fission products shown by black curves 
represent approximations based on experimental data by Prof. Sakamoto's group. This figure 
was done for us by Dr. Hiroshi Matsumura by adding our CEM03.01 results to Fig. 18 of the 
review |188j . 
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Figure 41: Excitation functions for the production of ^He and ^He from p+Fe. Left panel: 
Recent measurements by Ammon et al. |192j (filled squares) compared with our CEM03.03 
results (red lines), with results by TALYS [ 193j (blue line) from [192] . and with previous 
measurements shown with different symbols, as indicated (see references to old data in |192j ). 
Right panel: The same excitation functions as on the left, but predicted twelve years ago [76] 
with the CEM95 version of CEM (solid lines), compared with the experimental data available 
at that time (see references in [76] and [192] ). and with systematics by Konobeev and Korovin 
1 194] (thick long-dashed lines). Contributions from preequilibrium emission, evaporation, and 
from residual nuclei to the total CEM95 yields are shown by thin dashed lines, as indicated. 

indicated in Figs. 44 and 45. Fig. 44 shows the mass distribution of the product yields from 
the reaction 660 MeV p + ^^^I measured recently at JINR, Dubna [37] compared with our 
calculations with three recent versions of CEM. The standard version of CEM03.01 does not 
describe production of isotopes with mass number 26 < A < 63 from this reaction observed 
in the experiment [37]. These products are too heavy to be evaporated from compound nuclei 
and the target is too light to fission in GEM2, therefore not producing these isotopes as fission 
fragments (CEM03.01 and LAQGSM03.01 consider only "conventional" fission of preactinides 
and actinides and do not consider at all fission of nuclei with Z < 65). 

Fig. 45 shows a comparison [21] of the mass distributions of the yields of eight isotopes from 
Na to Mn produced in the reactions 1500, 1000, 750, 500, and 300 MeV/A ^'^Fe + p measured 
recently at GSI [162] with three recent versions of LAQGSM. We see that the standard version 
of LAQGSM03, referred as "03.01", fails to reproduce correctly production of fission-like heavy 
fragments from reactions with a medium-mass nuclear target at intermediate energies (see 
the solid lines on the bottom two panels of Fig. 45), just as do all previous versions, the 
standard CEM03.01 and its predecessors, and most other currently available models (see, e.g., 
|162t 1196] ). As already mentioned above, such nuclear targets are considered too light to 
fission in conventional codes (including GEM2 and all models currently employed in large-scale 
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transport codes). Similarly, the fragments are too light to be produced as spallation residues 
at these intermediate energies and too heavy to be produced via standard evaporation models. 
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Figure 42: Excitation functions for the production of ^^Ne and ^^Ar from p+Fe. Left panel: 
Recent measurements by Ammon et al. [192\ (filled squares) compared with our CEM03.03 
results (red lines), with results from TALYS |193j (blue solid line), and from INCL/ABLA 
[72l 1195] (dashed brown line) from |192j . and with previous measurements shown with different 
symbols, as indicated (see references to old data in [192] ). Right panel: The same excitation 
functions as on the left, but predicted twelve years ago [76] with the CEM95 version of CEM, 
compared with experimental data available at that time (see references in [76| and [192] ). as 
indicated. 

We note that the problem of a reliable description of intermediate-mass fragments from not 
too heavy targets is still unresolved not only in our standard versions of CEM and LAQGSM, 
but also in most other similar Monte-Carlo codes. Fig. 46 and 47 prove this with an example 
of comparison [191] the mass distributions of ^^Fe(p,x) reaction products measured at ITEP 
[191] and GSI |162l | at 300 and 1000 MeV energies with results from 15 different different 
code systems: MCNPX (INCL, CEM2k, BERTINI, ISABEL), LAHET (BERTINI, ISABEL), 
CEM03 (.01, .Gl, .SI), LAQGSM03 (.01, .Gl, .SI), CASCADE-2004, LAHETO, and BRIEFF 
(see references and details in |191] ). We see that all these models give a sufficiently good 
description of the mass yields of the products close to target nucleus mass {A > 35-40). In the 
mass range A < 30, however, a good description of the measured product nuclide yields is only 
given by the models that, apart from the conventional evaporation of nucleons and complex 
particles up to "^He, allow for evaporation of heavy clusters (the CEM and LAQGSM versions). 

We also note that none of these models gives a good quantitative description of the whole set 
of experimental data. However, as can be seen from Fig. 47, while still having these problems, 
our CEM and LAQGSM codes provide the lowest mean deviation factors between the calculated 
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and measured |162^ 1191] cross sections, averaged over all energies of the beams and detected 
products in comparison with other models (see more details in (1911). 




Tp (MeV) Tp (MeV) 

Figure 43: Excitation functions for the production of ^He, ^He, ^^Ne and ^^Ar from p+Ni. 
Recent measurements by Ammon et al. |192j (filled squares) compared with our CEM03.03 
results (red lines), with results from TALYS |193j (blue solid line), and from INCL/ABLA 
[72 | 1195] (dashed brown line) from |192j . and with previous measurements shown with different 
symbols, as indicated (see references to old data in |192j ). 

We have addressed this problem in two different ways [211 WJ\ 1196] : 

1) By implementing into CEM03.01 and LAQGSM03.01 the Statistical Multifragmentation 
Model (SMM) by Botvina et al. [155] . |197] - pUT] . to consider multifragmentation as a mode 
competitive to evaporation of particles and light fragments, when the excitation energy E* of 
a compound nucleus produced after the preequilibrium stage of a reaction is above 2x A MeV. 
This way, we have produced the "S" version of our codes ("S" stands for SMM), CEM03.S1 
and LAQGSM03.S1. 

CEM03.S1 and LAQGSM03.S1 are exactly the same as CEM03.01 and LAQGSM03.01, but 
consider also multifragmentation of excited nuclei produced after the preequilibrium stage of 
reactions, when their excitation energy is above 2 x A MeV, using SMM. In SMM, within 
the total accessible phase space, a micro-canonical ensemble of all breakup configurations, 
composed of nucleons and excited intermediate-mass fragments governs the disassembly of the 
hot remnant. The probability of different channels is proportional to their statistical weight. 
Several different breakup partitions of the system are possible. When after the preequilibrium 
stage of a reaction E* > 2 x A MeV and we "activate" SMM to calculate multifragmentation 
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in the 03. SI codes, the competitive evaporation process are calculated also with a version of 
the evaporation model by Botvina et al. from SMM, rather than using GEM2, as we do always 
in the standard 03.01 versions and in 03. SI when E* <2 'x A MeV and SMM is not invoked. 

A very detailed description of SMM together with many results obtained for many reactions 
may be found in |155j . |197] - [2(JT] and references therein; many details on SMM are presented 
in the first paper of Ref . [97] , therefore we do not repeat this here. 

2) By replacing the Generalized Evaporation Model GEM2 by Furihata [144j - [146j used in 
CEM03.01 and LAQGSM03.01 with the fission-like binary-decay model GEMINI of Charity 
|202] - [2U7] which considers production of all possible fragments. This way, we have produced 
the "G" version of our codes ("G" stands for GEMINI), CEM03.G1 and LAQGSM03.G1. 

CEM03.G1 and LAQGSM03.G1 are exactly the same as CEM03.01 and LAQGSM03.01, 
but use GEMINI instead of using GEM2. Within GEMINI, a special treatment based on the 
Hauser-Feshbach formalism is used to calculate emission of the lightest particles, from neutrons 
and protons up to beryllium isotopes. The formation of heavier nuclei than beryllium is mod- 
eled according to the transition-state formalism developed by Moretto [208] . All asymmetric 
divisions of the decaying compound nuclei are considered in the calculation of the probability 
of successive binary-decay configurations. GEMINI is described in details in the original publi- 
cations |202j ~ [208] : many details on GEMINI may be found also in the first paper of Ref. |97j . 
therefore we do not elaborate here. 
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Figure 44: Experimental [37] mass number distribution of the product yield from 660 MeV p 
+ 1291 compared with calculations by CEM03.01 (thick solid line), CEM03.S1 (thick dashed 
line), and CEM03.G1 (solid thin line), as indicated 
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Figure 45: Experimental [162] mass distributions of the yields of eight isotopes from Na to 
Mn produced in the reactions 1500, 1000, 750, 500, and 300 MeV/A ^^Fe + p compared with 
LAQGSM03 (solid lines on both panels), LAQGSM03+GEMINI (dashed lines, left panel), and 
LAQGSM03+SMM (dashed hues, right panel) results, respectively. 
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Figure 46: Mass distributions of ^^Fe(p,x) reaction products measured at ITEP (filled circles) 
^191j and GSI (open circles) |162] for 300 and 1000 MeV energies compared with results by 
14 codes shown with different lines as following: INCL/MCNPX (solid black), BRIEFF1.5.4g 
(dashed black), CEM03.01 (solid green), CEM2k/MCNPX (dashed green), CEM03.G1 (dotted 
green), BERTINI (MCNPX - solid blue, LAHET - dashed blue), ISABEL (MCNPX - sohd red, 
LAHET - dashed red, LAHETO - dotted red), LAQGSM03.01 (solid magenta), LAQGSM03.G1 
(dotted magenta), LAQGSM03.S1 (dashed-dotted magenta), CASCADE-2004 (cyan). This 
figure is adopted from Ref. [191] , where more details on this comparison may be found. 
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Figure 47: The predictive power of each tested code for the data shown in Fig. 46: the mean- 
square deviation factors between calculations and measured cross sections averaged over all 
energies and all products (see details in |191j ). 

Results by "S" and "G" versions of our codes are shown together with standard 03.01 results 
in Figs. 44, 45, and 46 (see also Figs. 28 and 29 in Section 6). We see that both "S" and "G" 
versions reproduce almost equally well the yields of intermediate-mass fragments from these 
reactions, that can not be described by the standard CEM03.01 and LAQGSM03.01. 

From one point of view, this is an achievement, as, with the "S" and "G" versions, we are 
able to describe reasonably well the reactions where the standard 03.01 versions do not work. On 
the other hand, this also makes the situation more intricate, as from these (and other similar) 
results we can not choose easily between the "S" and "G" versions, that makes more difficult to 
determine the real mechanisms of such reactions. We think that for such intermediate-energy 
proton-induced reactions the contribution of multifragmentation to the production of heavy 
fragments should not be very significant due to the relatively low excitation energies involved. 
Such fragments are more likely to be produced via the fission-like binary decays modeled by 
GEMINI. Multifragmentation processes are important and should be considered in reactions 
involving higher excitation energies; at excitations probably higher than the 2 MeV/nucleon 
considered here. We conclude [21] that it is impossible to make a correct choice between fission- 
like and fragmentation reaction mechanisms involved in these (or other similar) reactions merely 
by comparing model results with the measurements of only product cross sections; addressing 
this question will probably require analysis of two- or multi-particle correlation measurements. 
Our conclusion [2T] was confirmed by the recent coincidence measurements and calculations of 
residual and light particles in 1 GeV/A ^^Fe + p by Gentil et al. piMmH] . 

Several more examples of the difficult situation with determination of the mechanisms of 
fragment production from reactions induced by intermediate-energy projectiles using the "S" 
and "G" versions of our codes are presented below in Figs. 48 to 51. 

Fig. 48 shows data recently measured at GSI |211j on mass- and charge-product yield 
distributions and mean kinetic energies of all products as a function of the product mass 
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number from the reaction 1 GeV/nucleon ^^^Xe + p compared with calculations by the standard 
CEM03.03 and LAQGSM03.03, as well as by their "S" versions. More results for this reaction 
may be found in Ref. t212j. We see that the standard 03.03 versions of both GEM and LAQGSM 
codes fail to reproduce production of intermediate mass fragments with 20 < A < 70 from 
this reaction. With the "S" versions, we performed calculations for several different values of 
the excitations energy E* of the nuclei produced after the preequilibrium stage of reactions 
when we start to consider multifragmentation as a competitive to evaporation mechanism of 
nuclear reactions, namely, at E* > 2, 4, 4.5, and 5 MeV/nucleon in the case of GEM03.S1, 
and E* > 1.5, 2, and 4 MeV/nucleon for LAQGSM03.S1. We see that the best agreement of 
results by LAQGSM03.S1 (solid red lines in Fig. 48) with the GSI data |211j were achieved for 
E* > 2 MeV/nucleon, just as Dr. Botvina suggested to us and how it was implemented in both 
LAQGSM03.S1 and GEM03.S1 versions of our codes. But in the case of GEM03.S1 (sohd blue 
lines in Fig. 48), we got the best agreement with the data for a higher energy, namely E* > 4.5 
MeV/nucleon. 




Mass number A Charge number Z 




Figure 48: Mass- and charge-product yield distributions and mean kinetic energy of all products 
as functions of the product mass number from the reaction 1 GeV/nucleon ^^^Xe -|- p [211j (open 
circles) recently measured at GSI compared with calculations by the standard GEM03.03 and 
LAQGSM03.03 (dotted lines), as well as by their "S" versions (solid and dot-dashed lines) 
for different values of the excitations energy E* of the nuclei produced after the preequilibrium 
stage of reactions above which we start to consider multifragmentation as a competitive reaction 
mechanism, namely, at E* > 2, 4, 4.5, and 5 MeV/nucleon in the case of GEM03.S1, and 
E* > 1.5, 2, and 4 MeV/nucleon for LAQGSM03.S1, as indicated. 
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The fact that we get different "best" values of E* when we need to "activate" multifrag- 
mentation as a competition to the simple evaporation mechanism in different codes is not 
contradictory and does imply anything physical: The INC of CEM is completely different from 
the one of LAQGSM, therefore the mean mass < A > and charge number < Z > oi residual 
nuclei produced after the preequilibrium stage of a reaction with an excitation E* higher than 
a certain values, e.g., 2 MeV/nucleon, should be also different, as should be the distributions 
of such nuclei with respect to their E* . This results in turn in different fragments produced via 
multifragmentation from this reaction as predicted by CEM03.S1 and LAQGSM03.S1. If we 
would use another INC, e.g., INCL followed or not by preequilibrium emission, we would 
get still another "best" value of E* when we need to start considering multifragmentation, as 
another INC would predict other values of < A >, < Z >, and distributions of E* . This makes 
the situation of considering multifragmentation by different codes more intricate: the condition 
of its "activation" seems to be model-dependent. 

This situation becomes even more unclear if we refer to the results shown in Fig. 49; the 
same data [211] can be reproduced, even a little better, with the "G" versions of our codes, 
without considering multifragmentation at all. 
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Figure 49: The same experimental data [211] as in previous figure, but compared with results 
of calculations by CEM03.G1 (violet lines) and LAQGSM03.G1 (turquoise lines), as indicated. 



As described above, CEM03.G1 and LAQGSM03.G1 are exactly the same as CEM03.01 
and LAQGSM03.01, only replacing GEM2 [Ill]-[Tl6] with GEMINI [202]-[207]. As can be 
seen from Figs. 28, 29, 49, 50, and from results presented in Refs. [971 [HSl [212], the "G" 
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Figure 50: Comparison of all measured [213] cross sections of fission products from the reaction 
1 GeV/nucleon ^ospb on p (filled circles) with LAQGSM+GEMINI (LAQGSM03.G1) resuhs 
(open circles). Experimental data for isotopes of Ga and I are not available, so we present only 
our predictions. t_delay = 75 and sig_delay = 50 are used in GEMINI to calculate this reaction. 
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version of our codes allows us to describe well many fission and fragmentation reactions, es- 
pecially on targets below the actinide region, where the standard versions of our codes using 
GEM2 may not work well. 

We find that: 1) GEMINI merged with GEM or LAQGSM provides reasonably good results 
for medium-heavy targets without a fission delay time; 2) For preactinides, we have to use 
t_delay = 50-70 and sig^delay = 1-50, otherwise GEMINI provides too much fission — this may 
be related to the calculation of fission barriers of preactinides without strong ground-state shell 
corrections in the 2002 version of GEMINI we use; 3) The current version of GEMINI does not 
work well for low-energy fission of actinides (see Fig. 51). 




25 50 75 100 125 150 175 200 

Mass Number, A 

Figure 51: Mass distributions of fission fragments from the compound nucleus ^35^ ^ith an 
excitation of 14 MeV as predicted by GEMINI (solid line) and GEM2 (dashed line), respectively 
(adopted from [196]). No delay time in GEMINI for fission is used here. 

The physical basis of these shortcomings of GEMINI are understood and are discussed along 
with possible solutions to be employed in future versions of GEMINI by Dr. Gharity at this 
Workshop [Mi- 
ll. Summary 

Improved versions of the cascade-exciton model of nuclear reactions and of the Los Alamos 
quark-gluon string model have been developed recently at LANL and implemented in the codes 
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CEM03.01 and LAQGSM03.01, and in their slightly corrected 03.02 and 03.03 versions. The 
recent versions of our codes describe quite well and much better than their predecessors a large 
variety of nuclear reactions of interest to ADS, as well as reactions at higher energies of interest 
to NASA, accelerator shielding, astrophysics, and other applications, up to ~ 1 TeV/nucleon. 
(Energies above about 5 GeV are only accessible to LAQGSM.) 

What is more, our codes provide reasonable results even for low-energy reactions, where 
they are not easy to justify from a fundamental physics point of view (see e.g., [2Tj). 

We observe that codes developed for applications must not only describe reasonably well 
arbitrary reactions without any free parameters, but also not require too much computing 
time. Our own exercises on many different nuclear reactions with MCNPX [52] and LAHET 
|114[ 1131] transport codes using different event generators show that the current version of GEM 
is a little faster than the Bertini |173j and ISABEL |174] options for ING in these transport 
codes, and is several times faster than the INGL/ABLA j72 | 1195] option. The computing time 
of LAQGSM for heavy-ion reactions was tested recently by Dr. Gomez |214j with MGNPX 
in comparison with PHITS |215j . on simulations of a 400 MeV/nucleon uranium beam on a 
0.2-cm thick lithium target. He found that PHITS using the Quantum Molecular Dynamics 
(QMD) model to simulate heavy-ion reactions requires 210 times more computing time than 
LAQGSM does, making it very difficult if not impossible to use in applications requiring fast 
results, for example in real-time beam monitoring of positron emitting fragments used for PET 
tomography in heavy-ion cancer therapy; LAQGSM does not have these problems, providing 
quite fast and reliable results. 

Both our GEM and LAQGSM event generators still have some problems in a reliable de- 
scription of some fragments from intermediate-energy nuclear reactions on medium-mass nuclei, 
just as other similar modern Monte-Garlo codes do. To address these problems, we developed 
"S" and "G" modifications of our codes, that consider multifragmentation of nuclei formed after 
the preequilibrium stage of reactions when their excitation energy is above 2-5 x A MeV using 
the Statistical Multifragmentation Model (SMM) code by Botvina et al. |155] . |197] - [20T] . and 
the fission-like binary-decay model GEMINI by Gharity [202]-[207J, respectively. The "S" and 
"G" versions of our codes allow us to describe many fragmentation reactions that we are not 
able to reproduce properly with the standard version of our codes. However, there are still some 
problems to be solved in understanding the "true" mechanisms of fragment productions from 
many reactions, therefore we consider the present "S" and "G" versions of our codes only as 
working modifications and they are not yet implemented as event generators into our transport 
codes, in contrast to our "standard" 03.01, 03.02, and 03.03 versions used as event generators. 

The latest versions of our 03.01, 03.02, and 03.03 GEM and LAQGSM codes have been or 
are being incorporated as event generators into the transport codes MGNP6 |216] . MGNPX 
[S2], and MARS pT]. GEM03.01 was made available to the public via RSIGG at Oak Ridge and 
NEA/OEGD in Paris as the Gode Package PSR-0532. We also plan to make LAQGSM03.01 
and the latest 03.03 versions of our codes available to the public via RSIGG and NEA/OEGD 
in the future. 
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